WO2011121511A1 - Apparatus and method for forming a concentration image of the concentration of magnetic particles arranged in a field of view field of the invention - Google Patents

Abstract

The present invention relates to an apparatus and a method for forming a concentration image of the concentration of magnetic particles arranged in a field of view (28), in particular by applying the principles of magnetic particle imaging (MPI). Said apparatus comprises drive means, receiving means and reconstruction means (152) for reconstructing a concentration image of said field of view (28) from said acquired and sampled detection signals (s(t₁)) by the steps of: i) determining samples of a vector-valued function (Φ(t₁)) from said sampled detection signal (s(t₁)), and ii) reconstructing the position-dependent concentration (c*) of the magnetic particles in the field of view from the equation (I) wherein M is an integral operator modelling the data acquisition and p(t₁) are the known positions of the field-free point at the sampling times t₁ l = 0, 1, …, L.

Description

APPARATUS AND METHOD FOR FORMING A CONCENTRATION IMAGE OF THE CONCENTRATION OF MAGNETIC PARTICLES ARRANGED IN A FIELD OF VIEW FIELD OF THE INVENTION

The present invention relates to an apparatus and a method for forming a concentration image of the concentration of magnetic particles arranged in a field of view. Further, the present invention relates to a computer program for implementing said method on a computer and for controlling such an apparatus. The present invention relates particularly to the field of Magnetic Particle Imaging.

BACKGROUND OF THE INVENTION

Magnetic Particle Imaging (MPI) is an emerging medical imaging modality. The first versions of MPI were two-dimensional in that they produced two-dimensional images. Future versions will be three-dimensional (3D). A four-dimensional image of a non- static object can be created by combining a temporal sequence of 3D images to a movie, provided the object does not significantly change during the data acquisition for a single 3D image.

MPI is a reconstructive imaging method, like Computed Tomography or Magnetic Resonance Imaging. Accordingly, an MP image of an object's volume of interest is generated in two steps. The first step, referred to as data acquisition, is performed using an MPI scanner. The MPI scanner has means to generate a static magnetic gradient field, called the "selection field", which has a single field- free point (FFP) at the isocenter of the scanner. Moreover, this FFP is surrounded by a first sub-zone with a low magnetic field strength, which is in turn surrounded by a second sub-zone with a higher magnetic field strength. In addition, the scanner has means to generate a time-dependent, spatially nearly homogeneous magnetic field. Actually, this field is obtained by superposing a rapidly changing field with a small amplitude, called the "drive field", and a slowly varying field with a large amplitude, called the "focus field". By adding the time-dependent drive and focus fields to the static selection field, the FFP may be moved along a predetermined FFP trajectory throughout a

"volume of scanning" surrounding the isocenter. The scanner also has an arrangement of one or more, e.g. three, receive coils and can record any voltages induced in these coils. For the data acquisition, the object to be imaged is placed in the scanner such that the object's volume of interest is enclosed by the scanner's field of view, which is a subset of the volume of scanning.

The object must contain magnetic nanoparticles; if the object is an animal or a patient, a contrast agent containing such particles is administered to the animal or patient prior to the scan. During the data acquisition, the MPI scanner steers the FFP along a deliberately chosen trajectory that traces out the volume of scanning, or at least the field of view. The magnetic nanoparticles within the object experience a changing magnetic field and respond by changing their magnetization. The changing magnetization of the nanoparticles induces a time-dependent voltage in each of the receive coils. This voltage is sampled in a receiver associated with the receive coil. The samples output by the receivers are recorded and constitute the acquired data. The parameters that control the details of the data acquisition make up the "scan protocol".

In the second step of the image generation, referred to as image reconstruction, the image is computed, or reconstructed, from the data acquired in the first step. The image is a discrete 3D array of data that represents a sampled approximation to the position-dependent concentration of the magnetic nanoparticles in the field of view. The reconstruction is generally performed by a computer, which executes a suitable computer program. Computer and computer program realize a reconstruction algorithm. The reconstruction algorithm is based on a mathematical model of the data acquisition. As with all reconstructive imaging methods, this model can be formulated as an integral operator that acts on the acquired data; the reconstruction algorithm tries to undo, to the extent possible, the action of the model.

Such an MPI apparatus and method have the advantage that they can be used to examine arbitrary examination objects - e. g. human bodies - in a non-destructive manner and with a high spatial resolution, both close to the surface and remote from the surface of the examination object. Such an apparatus and method are generally known and have been first described in DE 101 51 778 Al and in Gleich, B. and Weizenecker, J. (2005),

"Tomographic imaging using the nonlinear response of magnetic particles" in Nature, vol. 435, pp. 1214-1217, in which also the reconstruction principle is generally described. The apparatus and method for magnetic particle imaging (MPI) described in that publication take advantage of the non-linear magnetization curve of small magnetic particles.

In conventional embodiments of MPI, the image is reconstructed by solving a generic, unstructured integral equation whose kernel is the so-called system function.

Reconstruction times tend to be long. Further, each time a new charge of magnetic particles is used, this system function has to be newly measured in a calibration procedure, which is very time consuming. In addition, the measured system function is degraded by measurement errors. SUMMARY OF THE INVENTION

It is an object of the present invention to provide an apparatus and a method for forming a concentration image of the concentration of magnetic particles arranged in a field of view, which enable a reconstruction of concentration images having a higher quality than conventionally reconstructed images, which enable a faster reconstruction, and which avoid the need to measure the system function in a calibration procedure.

In a first aspect of the present invention an apparatus for forming a

concentration image of the concentration of magnetic particles arranged in a field of view is presented comprising:

selection means comprising a selection field signal generator unit and selection field elements for generating a magnetic selection field having a pattern in space of its magnetic field strength such that a first sub-zone, having a low magnetic field strength and including a field- free point, and a second sub-zone having a higher magnetic field strength are formed in a field of view,

drive means comprising a drive field signal generator unit and drive field coils for changing the position in space of the two sub-zones in the field of view by means of a magnetic drive field so that the magnetization of the magnetic particles contained in said object changes locally,

receiving means comprising at least one signal receiving unit and at least one receiving coil for acquiring detection signals, which detection signals depend on the magnetization in the field of view, which magnetization is influenced by the change of the position in space of the first and second sub-zone, said detection signals being acquired in sampled form while the first sub-zone is moved along a trajectory,

reconstruction means for reconstructing a concentration image of said field of view from said acquired and sampled detection signals by the steps of:

i) determining samples of a vector- valued function Φ(ί) from said sampled detection signal, and

ii) reconstructing the position-dependent concentration c* of the magnetic particles in the field of view from the equation Φ(Π ) = (Me * ) (p(tt))— («ΑίΓ* )(ρ(ί₀)^'ϊ

wherein M is an integral operator modelling the data acquisition and p(ti) are the known positions of the field- free point at the sampling times ti 1 = 0, 1, ..., L.

In a further aspect of the present invention a corresponding method is presented, as well as a computer program for implementing said method.

Preferred embodiments of the invention are defined in the dependent claims. It shall be understood that the claimed method and the claimed computer program have similar and/or identical preferred embodiments as the claimed apparatus and as defined in the dependent claims.

The present invention provides a faster reconstruction algorithm. This algorithm solves a well structured integral equation of the convolution type and makes use of Non-uniform Fast Fourier Transforms. The acquired detection signals are preprocessed to make them satisfy this integral equation. The underlying MPI scanner should have certain features to make the reconstruction problem well solvable, as will be described below.

The reconstruction means can generally be implemented in various ways, e.g. by hard- and/or software, for instance by use of a processor.

The integral operator M represents an indirect model of the data acquisition, which is derived from a direct model by deconvolution of the impulse response of the receiving means (which deconvolution converts detection signals to electrical voltages), integration with respect to time (which converts electrical voltages to magnetic fluxes), and replacement of the time variable by a position variable (which almost generates a spatial convolution operator). These steps will be explained in more detail below.

Hence, in preferred embodiments said reconstruction means is adapted for reconstructing the samples of said vector-valued function from said detection signals by discrete deconvolution and integration, and for reconstructing the concentration image from the said samples of the said deconvolved and integrated detection signals, preferably by use of a conjugate gradient method.

For the definition of the integral operator M various options exist that allow a faster reconstruction of the particle concentration.

In one embodiment the integral operator M is defined by V wherein V denotes the field of view, R(x) is the sensitivity matrix of the at least one receiving coil, H*(x, y) is the total magnetic field at position x when the field-free point is at position y, and P(h) is a vector valued function defined by

P(h) := ο /ΗΙΙ βίΗ),

wherein μο is the magnetic permeability of vacuum, P(H) is the magnetization curve, and the vector e(h) is the unit vector along the direction of the vector h. Thereby,^— J is the vector whose components are the magnetic fluxes through the receive coils caused by the concentration f of magnetic particles when the field-free point is at position y.

More advantageous embodiments are obtained by exploiting the hidden convolutional structure of the operator M . Indeed, using the approximation

H ^' *? }') ^:¾ (x) — · H ^: ) _{s w}herein H^s(x) is the selection field at position x, converts the operator M to the intermediate operator

(M₂f)(y) :^ I R{x)^rP(H^s(y) - H^s (x})/(x) < x , y€ W.

J v

where W denotes the volume of scanning. The operator M ₂ may now be easily transformed to an operator with a convolutional structure. The required transformation is given by the mapping x i— > H^s(x) . The transformed operator can be written in the form

(¾/)(s) :== / R(r)^rP(s - r)/(r) d i s e H^S ( W), where ^ ^*j ^( (H ) (*^*) ) is the transformed sensitivity matrix. The transformed function f(r) is related to the function f(x) by f{r) : J(_Hs^ (r) (H^S )^_1 (r) , where ^J (H^s j^{"1 iSt the J acobian determinant of the} inverse mapping

A reconstruction algorithm designed for 3 would produce an image expressed in the transformed coordinate r, which ranges in the space H^s(v). Therefore, one would have to transform the reconstructed image back to spatial coordinates.

Alternatively, in a further embodiment the operator M is replaced by the operator Ji₃ defined by

(«& _¾/)< j-) : = (x)^TY(H^&(y - x))f(x) d x_f

This embodiment requires that the magnetic selection field be substantially linear (so that H^s (y) - H^s (x) = H^s (y - x) at least approximately), but it does not involve a coordinate transformation.

The operators =·^ 3 and Jt₃ are advantageous in that their convolutional structure allows for faster reconstruction algorithms that make use of Non-uniform Fast Fourier Transforms.

In a further embodiment, the operator --^ 3 is replaced by the operator <^ defined by

(¾ )(») ^:= Rj I P{$ - r)/(r) i wherein Ro is a constant sensitivity matrix. This embodiment requires that this sensitivity matrix be substantially independent of x, i.e. that the magnetic selection field is substantially linear.

If, in addition, the selection field is substantially linear, the operator *¾4 may be further replaced by the operator Jt₄ defined by

f

" ¹ '·· ^»·

where Ro is again a constant sensitivity matrix.

Using the operators <^4 and Ji₄ allows even faster reconstruction algorithms, since these operators and Ji₄ are pure convolution operators, whereas the operators -^ i and Ji₃ are generalized convolution operators.

In addition to the selection means and the drive means for generating magnetic fields, focus means are provided in an embodiment for changing the position in space of the field- free point by means of a magnetic focus field. Such a focus field should have the same (or a similar) spatial distribution as the drive field. Separate or the same (preferably) means (e.g. coils) can be used as the focus means and the drive means. If the same means (in particular coils) are used for generating the focus field and the drive field, these fields have the same spatial distribution. Equal or at least similar spatial distributions of the drive and focus fields are required to be able to relate time t with the position of the field- free point at that time, p(t). If the spatial distributions of these fields are very different, the method may no longer be useful or lead to unfavorable results. The basic difference between the drive field and the focus field is that the frequencies are much lower (e.g. <1 kHz, typically <100 Hz) for the focus field than for the drive field, but the amplitudes of the focus field are much higher (e.g. 200 mT compared to 20 mT for the drive field). These fields are used to move the FFP to a desired position. The drive field is required in addition to the focus field since the detection signal obtainable with only the focus field would not be usable for the desired purpose, as the frequencies produced in the object of interest are much too low (typically <10 kHz).

BRIEF DESCRIPTION OF THE DRAWINGS

These and other aspects of the invention will be apparent from and elucidated with reference to the embodiment(s) described hereinafter. In the following drawings

Fig. 1 shows a first embodiment of an MPI apparatus,

Fig. 2 shows an example of the selection field pattern produced by an apparatus as shown in Fig. 1,

Fig. 3 shows a second embodiment of an MPI apparatus,

Fig. 4A shows a block diagram of an embodiment of an MPI apparatus according to the present invention,

Fig. 4B shows a block diagram of an embodiment of an MPI apparatus, Fig. 5 shows an FFP trajectory for the first embodiment of the MPI apparatus,

Fig. 6 shows an example of an FFP trajectory to be created by the focus field, and

shows a graph of the Langevin function.

DETAILED DESCRIPTION OF THE INVENTION

Before the details of the present invention shall be explained, basics of magnetic particle imaging shall be explained in detail with reference to Figs. 1 to 4. In particular, two embodiments of an MPI scanner for medical diagnostics will be described. An informal description of the data acquisition will also be given. The similarities and differences between the two embodiments will be pointed out.

The first embodiment 10 of an MPI scanner shown in Fig. 1 has three pairs 12, 14, 16 of coaxial parallel circular coils, these coil pairs being arranged as illustrated in Fig. 1. These coil pairs 12, 14, 16 serve to generate the selection field as well as the drive and focus fields. The axes 18, 20, 22 of the three coil pairs 12, 14, 16 are mutually orthogonal and meet in a single point, designated the isocenter 24 of the MPI scanner 10. In addition, these axes 18, 20, 22 serve as the axes of a 3D Cartesian x-y-z coordinate system attached to the isocenter 24. The vertical axis 20 is nominated the y-axis, so that the x- and z-axes are horizontal. The coil pairs 12, 14, 16 are named after their axes. For example, the y-coil pair 14 is formed by the coils at the top and the bottom of the scanner. Moreover, the coil with the positive (negative) y-coordinate is called the y -coil (y -coil), and similarly for the remaining coils. When more convenient, the coordinate axes and the coils shall be labelled with x_{l s} x₂, and x₃, rather than with x, y, and z.

The scanner 10 can be set to direct a predetermined, time-dependent electric current through each of these coils 12, 14, 16, and in either direction. If the current flows clockwise around a coil when seen along this coil's axis, it will be taken as positive, otherwise as negative. To generate the static selection field, a constant positive current I^s is made to flow through the z⁺-coil, and the current -I^s is made to flow through the z -coil. The z-coil pair 16 then acts as an anti-parallel circular coil pair.

The magnetic selection field, which is generally a magnetic gradient field, is represented in Fig. 2 by the field lines 50. It has a substantially constant gradient in the direction of the (e.g. horizontal) z-axis 22 of the z-coil pair 16 generating the selection field and reaches the value zero in the isocenter 24 on this axis 22. Starting from this field- free point (not individually shown in Fig. 2), the field strength of the magnetic selection field 50 increases in all three spatial directions as the distance increases from the field-free point. In a first sub-zone or region 52 which is denoted by a dashed line around the isocenter 24 the field strength is so small that the magnetization of particles present in that first sub-zone 52 is not saturated, whereas the magnetization of particles present in a second sub-zone 54 (outside the region 52) is in a state of saturation. In the second sub-zone 54 (i.e. in the residual part of the scanner's field of view 28 outside of the first sub-zone 52) the magnetic field strength of the selection field is sufficiently strong to keep the magnetic particles in a state of saturation.

By changing the position of the two sub-zones 52, 54 (including the field- free point) within the field of view 28 the (overall) magnetization in the field of view 28 changes. By determining the magnetization in the field of view 28 or physical parameters influenced by the magnetization, information about the spatial distribution of the magnetic particles in the field of view 28 can be obtained. In order to change the relative spatial position of the two sub-zones 52, 54 (including the field-free point) in the field of view 28, further magnetic fields, i.e. the magnetic drive field, and, if applicable, the magnetic focus field, are superposed to the selection field 50.

To generate the drive field, a time dependent current I°i is made to flow through both x-coils 12, a time dependent current I^D ₂ through both y-coils 14, and a time dependent current I^D ₃ through both z-coils 16. Thus, each of the three coil pairs acts as a parallel circular coil pair. Similarly, to generate the focus field, a time dependent current I^Fi is made to flow through both x-coils 12, a current I^F ₂ through both y-coils 14, and a current I^F ₃ through both z-coils 16.

It should be noted that the z-coil pair 16 is special: It generates not only its share of the drive and focus fields, but also the selection field. The current flowing through the z^±-coil is I^D ₃ + I^F ₃ ± I^s. The current flowing through the remaining two coil pairs 12, 14 is I°_k + I^F _k, k = 1, 2. Because of their geometry and symmetry, the three coil pairs 12, 14, 16 are well decoupled. This is wanted.

Being generated by an anti-parallel circular coil pair, the selection field is rotationally symmetric about the z-axis, and its z-component is nearly linear in z and independent of x and y in a sizeable volume around the isocenter 24. In particular, the selection field has a single field- free point (FFP) at the isocenter. In contrast, the

contributions to the drive and focus fields, which are generated by parallel circular coil pairs, are spatially nearly homogeneous in a sizeable volume around the isocenter 24 and parallel to the axis of the respective coil pair. The drive and focus fields jointly generated by all three parallel circular coil pairs are spatially nearly homogeneous and can be given any direction and strength, up to some maximum strength. The drive and focus fields are also time- dependent. The difference between the focus field and the drive field is that the focus field varies slowly in time and may have a large amplitude, while the drive field varies rapidly and has a small amplitude. There are physical and biomedical reasons to treat these fields differently. A rapidly varying field with a large amplitude would be difficult to generate and hazardous to a patient.

The embodiment 10 of the MPI scanner has at least one further pair, preferably three further pairs, of parallel circular coils, again oriented along the x-, y-, and z- axes. These coil pairs, which are not shown in Fig. 1, serve as receive coils. As with the coil pairs 12, 14, 16 for the drive and focus fields, the magnetic field generated by a constant current flowing through one of these receive coil pairs is spatially nearly homogeneous within the field of view and parallel to the axis of the respective coil pair. The receive coils are supposed to be well decoupled. The time-dependent voltage induced in a receive coil is amplified and sampled by a receiver attached to this coil. More precisely, to cope with the enormous dynamic range of this signal, the receiver samples the difference between the received signal and a reference signal. The transfer function of the receiver is non-zero from zero Hertz ("DC") up to the frequency where the expected signal level drops below the noise level.

The embodiment 10 of the MPI scanner shown in Fig. 1 has a cylindrical bore 26 along the z-axis 22, i.e. along the axis of the selection field. All coils are placed outside this bore 26. For the data acquisition, the patient (or object) to be imaged is placed in the bore 26 such that the patient's volume of interest - that volume of the patient (or object) that shall be imaged - is enclosed by the scanner's field of view 28 - that volume of the scanner whose contents the scanner can image. The patient (or object) is, for instance, placed on a patient table. The field of view 28 is a geometrically simple, isocentric volume in the interior of the bore 26, such as a cube, a ball, or a cylinder. A cubical field of view 28 is illustrated in Fig. 1.

The size of the first sub-zone 52 is dependent on the strength of the gradient of the magnetic selection field and on the field strength of the magnetic field required for saturation, which in turn depends on the magnetic particles. For a sufficient saturation of typical magnetic particles at a magnetic field strength of 80 A/m and a gradient (in a given space direction) of the field strength of the magnetic selection field amounting to 50xl0³ A/m², the first sub-zone 52 in which the magnetization of the particles is not saturated has dimensions of about 1 mm (in the given space direction).

The patient's volume of interest is supposed to contain magnetic nanoparticles. Prior to the diagnostic imaging of, for example, a tumor, the magnetic particles are brought to the volume of interest, e.g. by means of a liquid comprising the magnetic particles which is injected into the body of the patient (object) or otherwise administered, e.g. orally, to the patient.

An embodiment of magnetic particles comprises, for example, a spherical substrate, for example, of glass which is provided with a soft-magnetic layer which has a thickness of, for example, 5 nm and consists, for example, of an iron-nickel alloy (for example, Permalloy). This layer may be covered, for example, by means of a coating layer which protects the particle against chemically and/or physically aggressive environments, e.g. acids. The magnetic field strength of the magnetic selection field 50 required for the saturation of the magnetization of such particles is dependent on various parameters, e.g. the diameter of the particles, the used magnetic material for the magnetic layer and other parameters. In the case of e.g. a diameter of 10 μιη, a magnetic field of approximately 800 A/m (corresponding approximately to a flux density of 1 mT) is then required, whereas in the case of a diameter of 100 μιη a magnetic field of 80 A/m suffices. Even smaller values are obtained when a coating of a material having a lower saturation magnetization is chosen or when the thickness of the layer is reduced. Magnetic particles that can generally be used are commercially available under the trade name Resovist.

For further details of the generally usable magnetic particles and particle compositions, the corresponding parts of EP 1304542, WO 2004/091386, WO 2004/091390, WO 2004/091394, WO 2004/091395, WO 2004/091396, WO 2004/091397, WO 2004/091398, WO 2004/091408 are herewith referred to, which are herein incorporated by reference. In these documents more details of the MPI method in general can be found as well.

During the data acquisition, the x-, y-, and z-coil pairs 12, 14, 16 generate a position- and time-dependent magnetic field, the applied field. This is achieved by directing suitable currents through the field generating coils. In effect, the drive and focus fields push the selection field around such that the FFP moves along a preselected FFP trajectory that traces out the volume of scanning - a superset of the field of view. The applied field orientates the magnetic nanoparticles in the patient. As the applied field changes, the resulting magnetization changes too, though it responds nonlinearly to the applied field. The sum of the changing applied field and the changing magnetization induces a time-dependent voltage V_k across the terminals of the receive coil pair along the Xk-axis. The associated receiver converts this voltage to a signal S_k, which it processes further.

Like the first embodiment 10 shown in Fig. 1 , the second embodiment 30 of the MPI scanner shown in Fig. 3 has three circular and mutually orthogonal coil pairs 32, 34, 36, but these coil pairs 32, 34, 36 generate the selection field and the focus field only. The z- coils 36, which again generate the selection field, are filled with ferromagnetic material 37. The z-axis 42 of this embodiment 30 is oriented vertically, while the x- and y-axes 38, 40 are oriented horizontally. The bore 46 of the scanner is parallel to the x-axis 38 and, thus, perpendicular to the axis 42 of the selection field. The drive field is generated by a solenoid (not shown) along the x-axis 38 and by pairs of saddle coils (not shown) along the two remaining axes 40, 42. These coils are wound around a tube which forms the bore. The drive field coils also serve as receive coils.

To give a few typical parameters of such an embodiment: The z-gradient of the selection field, G, has a strength of G/μο = 2.5 T/m, where μο is the vacuum permeability. The temporal frequency spectrum of the drive field is concentrated in a narrow band around 25 kHz (up to approximately 100 kHz). The useful frequency spectrum of the received signals lies between 50 kHz and 1 MHz (eventually up to approximately 10 MHz). The bore has a diameter of 120 mm. The biggest cube 28 that fits into the bore 46 has an edge length of 120 mm/ ~ 84 mm.

Since the construction of field generating coils is generally known in the art, e.g. from the field of magnetic resonance imaging, this subject need not be further elaborated herein.

In an alternative embodiment for the generation of the selection field, permanent magnets (not shown) can be used. In the space between two poles of such

(opposing) permanent magnets (not shown) there is formed a magnetic field which is similar to that shown in Fig. 2, that is, when the opposing poles have the same polarity. In another alternative embodiment, the selection field can be generated by a mixture of at least one permanent magnet and at least one coil.

Fig. 4A shows a general block diagram of an MPI apparatus 100 according to the present invention. The general principles of magnetic particle imaging explained above are valid and applicable to this embodiment as well, unless otherwise specified.

The embodiment of the apparatus 100 shown in Fig. 4A comprises various sets of coils for generating the desired magnetic fields. First, the coils and their functions in MPI shall be explained.

For generating the magnetic selection field explained above, selection means are provided comprising a set of selection field coils 116, preferably comprising at least one pair of coil elements. The selection means further comprises a selection field signal generator unit 110. Preferably, a separate generator subunit is provided for each coil element (or each pair of coil elements) of the set 116 of selection field coils. Said selection field signal generator unit 110 comprises a controllable selection field current source 112 (generally including an amplifier) and a filter unit 114 which provide the respective section field coil element with the selection field current to individually set the gradient strength of the selection field. However, since the selection field is generally static, the filter unit 114 can also be omitted. Preferably, a constant current is provided. If the selection field coil elements are arranged as opposite coils, e.g. on opposite sides of the field of view, the selection field currents of the opposite coils are preferably oppositely oriented.

The selection field signal generator unit 110 can be controlled by a control unit 150, which preferably controls the selection field current generation 110 such that the sum of the field strength and the sum of the gradient strength of all spatial fractions of the selection field is maintained at a predefined level. For this purpose the control unit 150 can also be provided with control instructions by a user according to the desired application of the MPI apparatus, which, however, is preferably omitted according to the present invention.

For the generation of a magnetic focus field the apparatus 100 further comprises focus means comprising a set of focus field coils, preferably comprising three pairs 126a, 126b, 126c of oppositely arranged focus field coil elements. Said magnetic focus field is generally used for changing the position in space of the first and second sub-zones. The focus field coils are controlled by a focus field signal generator unit 120, preferably comprising a separate focus field signal generation subunit for each coil element (or at least each pair of coil elements) of said set of focus field coils. Said focus field signal generator unit 120 comprises a focus field current source 122 (preferably comprising a current amplifier) and a filter unit 124 for providing a focus field current to the respective coil of said subset of coils 126a, 126b, 126c which shall be used for generating the magnetic focus field. The focus field current unit 120 is also controlled by the control unit 150. With the present invention, the filter unit 124 may also be omitted.

For generating the magnetic drive field the apparatus 100 further comprises drive means comprising a subset of drive field coils, preferably comprising three pairs 136a, 136b, 136c of oppositely arranged drive field coil elements. The drive field coils are controlled by a drive field signal generator unit 130, preferably comprising a separate drive field signal generation subunit for each coil element (or at least each pair of coil elements) of said set of drive field coils. Said drive field signal generator unit 130 comprises a drive field current source 132 (preferably including a current amplifier) and a filter unit 134 (which may also be omitted with the present invention) for providing a drive field current to the respective drive field coil. The drive field current source 132 is adapted for generating a time- dependent current and is also controlled by the control unit 150.

It should be noted that in the embodiment of the apparatus 10 shown in Fig. 1 identical coils are preferably used for generating the magnetic drive field and the magnetic focus field.

For the detection of the signals, detection receiving means 148, in particular a receiving coil, and a signal receiving unit 140, which receives signals detected by said receiving means 148, are provided. Preferably, three receiving coils 148 and three receiving units 140 - one per receiving coil - are provided in practice, but more than three receiving coils and receiving units can be also used, in which case the acquired detection signals are not 3-dimensional but K-dimensional, with K being the number of receiving coils. According to the embodiment shown in Fig. 4A, a differential signal D_k is obtained by a subtraction unit 142, which subtracts approximation detection signals S_k from the actual detection signals S_k obtained by the receiving coil 148 connected to the k-th receiving unit 140. The approximation detection signals S_k are generally stored in digital form in an approximation detection signal storage unit 160. Details of how the

approximation detection signals S_k are obtained, will be explained below. The obtained differential signal D_k is then converted into a digital signal by the analog-to-digital converter 144 and stored in a storage unit 162. Alternatively, the differential signal D_k may be stored directly in a reconstruction unit 152.

In the same way, but at a different time, another differential signal D_k° is obtained as the difference between the approximation detection signal S_k and a calibration detection signal S_k°. The calibration detection signal S_k° is obtained by the signal receiving coil 148 with an empty scanner, i.e. without magnetic particles contained in the field of view. This differential signal D_k° is also converted into a digital signal by the analog-to-digital converter 144 and stored in the storage unit 162. Alternatively, the differential signal D_k° may be stored directly in the reconstruction unit 152.

It shall be noted that for the generation and sampling of the differential signals D_k and D_k° separate subtraction units and separate analog-to-digital converters can be used, differently from the embodiment shown in Fig. 4A.

Finally, the two differential signals D_k and D_k° are subtracted from each other in the reconstruction unit 152, thus resulting in the corrected detection signals Sk.

Alternatively, a separate subtraction unit for subtracting the differential signal D_k° from the differential signal D_k may be provided.

These corrected detection signals are then used by the reconstruction unit (also called reconstruction means) 152, which reconstructs the spatial distribution of the magnetic particles from these signals and the trajectory along which the field- free point moved during the data acquisition and which the reconstruction unit 152 obtains from the control unit 150. The reconstructed spatial distribution of the magnetic particles is finally transmitted via the control means 150 to a computer 154, which displays it on a monitor 156. Thus, an image can be displayed showing the distribution of magnetic particles in the field of view of the MPI scanner.

Further, an input unit 158 may be provided, for example a keyboard. A user is therefore able to choose the scan protocol and to set other parameters. According to the embodiment shown in Fig. 4A the reconstruction unit 152 (i.e. the reconstruction means) for reconstructing a concentration image of said field of view from corrected detection signals is adapted for determining said corrected detection signals from i) actual detection signals acquired at different positions of the first sub-zone with the magnetic particles arranged in the field of view, ii) calibration detection signals acquired at different positions of the first sub-zone with no magnetic particles arranged in the field of view, and iii) approximation detection signals determined as an approximation to said calibration detection signals, wherein said corrected detection signals are obtained by subtracting the difference between said calibration detection signals and said approximation detection signals from the difference between said actual detection signals and said approximation detection signals.

The applied total magnetic field Ho orientates the magnetic nanoparticles in the field of view. This results in a magnetization M and a slightly modified total magnetic field H, which together make up the magnetic flux density B = μο (H + M). The temporal change of the magnetic flux density induces a time dependent voltage V_k in the receive coil pair(s) 148 along the Xk-axis, k = 1, 2, 3. These voltages may be written as

Vk = V_k° + Vk, k = l, 2, 3,

where V_k° and Vk are the voltages induced by H and by M, respectively. For the

reconstruction, it is required to relate the voltages Vk induced by M with the magnetic nanoparticles causing M.

The magnetization-induced voltages are related to the acquired signals. This relationship can be written as

S_k = at (b_k · V_k) = at (b_k · (Vk° + v_k)), k = 1 , 2, 3,

where the time-dependent and known function bk is the impulse response of the receiver associated with the receive coil along the Xk-axis, and the unknown constant <¾ is proportional to the gain of this receiver. Assuming properly adjusted gains, it may be assumed that the scanner (i.e. the apparatus) acquires the signals (also called "actual detection signals")

S_k = a(bk- (Vk° + Vk)), k = l, 2, 3,

where a is a single unknown constant. With some extra effort it is possible to determine a, too, but it can be left unknown. For the model and the reconstruction, the signals

Sk = a (b_k-Vk), k = l, 2, 3, are desired. These signals, however, cannot be measured directly. However, by repeating the data acquisition with an empty MPI scanner, i.e. without any magnetic particles and an object in the field of view, the signals (also called "calibration detection signals")

S_k° = a (b_k · V_k°), k = l, 2, 3,

can be acquired. Strictly speaking, the signals thus obtained are induced not by the total field H, but by the applied field H₀. However, the difference between H and H₀ is negligible. The wanted signal Sk is then found by subtracting S_k° from S_k,

Sk = Sk - Sk°, k = l, 2, 3.

One problem remains: The signals S_k and S_k° have an enormous dynamic range and ordinary A/D converters cannot sample these signals with the required accuracy. To solve this problem, a crude approximation S_k (also called "approximation detection signals") to S_k° is obtained and the difference signals ("intermediate differential signals") D_k

* 0 0 *

= S_k - S_k and D_k = S_k - S_k are sampled, which have a lower dynamic range and hence can be sampled using ordinary A/D converters. The wanted signal ("corrected detection signal") is the difference of the difference signals:

s_k = (S_k - S_k ^*) - (Sk° - S_k ^*), k = l, 2, 3.

Suitable signals S_k can be predetermined by measurement or by simulation. In the case of simulation, a digital version of S_k is obtained by a separate (external) processing unit or by the reconstruction unit 152. In such a simulation, the Biot-Savart law describing the magnetic field generated by an electric current may be used.

In case of measurement, a digital version of S_k is obtained by m-bit AD conversion of the analog calibration detection signals S_k°, e.g. by use of an m-bit AD converter. This can be also be done by the apparatus shown in Fig. 4A, namely by measuring signals without any object (patient) being placed in the field of view and having the subtraction unit 142 subtract nothing from the measured signals.

A block diagram of an alternative embodiment of an MPI apparatus 100 is shown in Fig. 4B. Most elements of said embodiment are identical to elements used in the embodiment shown in Fig, 4A, are identified by identical numerals and will not be explained again. Mainly the differences to the embodiment shown in Fig. 4A will be explained.

According to this embodiment the signal receiving unit 140 comprises a filter unit 143 for filtering the received detection signals. The aim of this filtering is to separate measured values, which are caused by the magnetization in the examination area which is influenced by the change in position of the two part-regions 52, 54, from other, interfering signals. To this end, the filter unit 143 may be designed for example such that signals which have temporal frequencies that are smaller than the temporal frequencies with which the receiving coil 148 is operated, or smaller than twice these temporal frequencies, do not pass the filter unit 1. The signals are then transmitted via an amplifier unit 142 to an analog/digital converter 144 (ADC). The digitalized signals produced by the analog/digital converter 144 are fed to the reconstruction unit 152. However, the filter also suppresses valuable information about the image to be reconstructed. In contrast, the embodiment according to Fig. 4 A does not suppress such information.

For the following description, an apparatus 100 according to Fig. 4A is considered, which represents a possible embodiment of an MPI apparatus according to the present invention. The same idea described in the following may, however, also be applied in the embodiment shown in Fig. 4B, both with or without (preferred) the filter unit 143. The reconstruction unit 152 is adapted according to the present invention for reconstructing a concentration image of the field of view 28 from the acquired samples s(ti), 1 = 0, 1, ..., L, of the vector- valued detection signal s(t) = (si(t), s₂(t), s₃(t)). This is done in two steps. In a first step, samples Φ(ίι), 1 = 1, ..., L, of a vector- valued function < (t) are determined from said samples s(ti) of the detection signal s(t). Thereafter, in a second step the position-dependent concentration c* of the magnetic particles in the field of view is reconstructed from a discrete version of the equation

wherein M is an integral operator modelling the data acquisition and p(ti) is the position of the field- free point at sampling time t_ls 1 = 0, 1, ..., L.

The operator M has the property that the k-th component of ^— H}^') represents the magnetic flux through the receive coil along the Xk-axis caused by the concentration f(x) of magnetic particles when the field- free point is at position y. The operator M is related to the vector- valued function < (t) by the above equation. Accordingly the vector-valued function < (t) is proportional to the difference between the magnetic flux at time t (or t_ls if samples are considered rather than continuous functions) and the magnetic flux at time to.

As pointed out earlier, the acquired detection signal Sk(t) is the convolution of the voltage Vk(t) and the impulse response b_k(t) of the receiver associated with the receiving coil along the Xk-axis. This impulse response is deconvolved from the voltage by convolving the signal with the inverse filter b_k*(t) The inverse filter exists since the amplifiers of the receiving coils have a transfer function without zeros in the range from zero Hertz up to an appropriate maximum upper frequency. Thereafter, the voltage Vk(t) is integrated with respect to time. The result is a magnetic flux, up to an additive integration constant.

The concentration image may be reconstructed by discretizing the above equation and solving the resulting system of equations using, for example, the conjugate gradient method. However, because the system matrix is huge, dense, and not well structured, the computational cost would be high. For the same reasons, alternative computational methods would also be relatively slow. But as pointed out earlier, the operator M can be replaced by various operators that have a convolutional structure. Using one of these operators gives the system matrix a special structure which makes it possible to speed up the conjugate gradient method using Non-uniform Fast Fourier Transforms.

In the following a more detailed mathematical description of essential aspects of the present invention is provided.

In the following, an explicit mathematical model of the data acquisition of Magnetic Particle Imaging (MPI) is derived. This model is an integral operator whose kernel can be expressed in terms of the various magnetic fields generated by the MPI scanner, the magnetization curve of the magnetic particles, the sensitivity patterns of the receive coils, and the impulse responses of the receivers. Once these ingredients have been measured or determined otherwise, the kernel is known. Moreover, the model is closely related to a convolution operator. This fact makes it possible to discuss the solvability of the reconstruction problem and also paves the way for fast reconstruction algorithms based on uniform and nonuniform fast Fourier transforms. Several such reconstruction algorithms are outlined. To make the reconstruction problem well solvable, the MPI scanner should have certain physical/technical properties. These properties are pointed out.

List of Abbreviations

CT Computed Tomography

DC Direct Current

FFP Field Free Point

FFT Fast Fourier Transform

iCD Clinical Demonstrator; nickname of the example MPI scanner iFFT inverse Fast Fourier Transform

MPI Magnetic Particle Imaging

MRI Magnetic Resonance Imaging

NFFT Nonuniform Fast Fourier Transform

PCD Preclinical Demonstrator

SI Systeme International (d' Unites) List of Symbols

R the set of real numbers

time

r_s, r_e start and end times of the data acquisition

X, y points in three-dimensional space

Xi , x-2. A' 3 Cartesian components of x

x% ^*, z ditto

/½ magnetic permeability of vacuum

c(x) concentra tion of magnetic nanoparticles

p( ) FFP trajectory

s(f) vector formed by the received signals

( deconvolved and i ntegrated signal vector

F field of view, contains magnetic nanoparticles

W volume of scanning, traced by the FFP trajectory

I.^D( vector formed by the currents flowing through the drive held coils l^F(t) vector formed by the currents flowing through the focus field coils

1(f) sum of I^D(?) and l^F(f); moves the FFP to p(/)

G(x. t) system function

C(x) field matrix of the held generating coils

R(x) sensitivity matrix of the receive coils

H^D(x, r ) drive field

H^F(x, f ) focus field

H^s(x) selection field

H(x. t) total magnetic held, sum of selection, drive, and focus fields e(H) unit: vector along the direction of H

I*(x) current map, defined by I*(x) — C(x)^{~ 1}H^s(x)

H* (x, ) the magnetic held defined by H^*(x, y) H^s(x) + C(x)I* (y) p(H ) magnetization curve

P(H) vector-valued variant of P(H\ defined by P(H) /¼.P(|[H.jj)e(H

M(x, ) magnetization

M, Mi , ... integral operators that model the data acquisition

Kfx, y) kernel of M 1 Introduction

Magnetic Particle Imaging (MPI) is an emerging medical imaging modality under development at Philips Research Europe. The first versions of MPI were two- dimensional in that they produced two-dimensional images. Future versions, including the Preclinical Demonstrator (PCD) under construction at Philips Research Europe, will be three-dimensional (3D). This note concentrates on 3D MPI of static objects. A time-dependent, or 4D, image of a non-static object can be created by combining a temporal sequence of 3D images to a movie, provided the object does not significantly change during the acquisition of the data for a single 3D image.

MPI is a reconstructive imaging method, like Computed Tomography (CT) or Magnetic Resonance Imaging (MRI). Accordingly, an MP image of an object's volume of interest is generated in two steps. The first step, referred to as data acquisi- tion, is performed using an MPI scanner. The MPI scanner has means to generate a static magnetic gradient field, the selection field, which has a single field free point (FFP) at the isocenter of the scanner. In addition, the scanner has means to generate a time-dependent, spatially nearly homogeneous magnetic field. Actually, this field is obtained by superposing a rapidly changing field with a small amplitude, the drive field, and a slowly varying field with a large amplitude, the focus field. By adding the time-dependent drive and focus fields to the static selection field, the FFP may be moved along a predetermined FFP trajectory throughout a volume of scanning around the isocenter. The scanner also has an arrangement of three receive coils and can record any voltages induced in these coils. For the data acquisition, the object to be imaged is placed in the scanner such that the object's volume of interest is enclosed by the scanner's field of view, which is a subset of the volume of scanning. The object must contain magnetic nanoparticles; if the object is an animal or a patient, a contrast agent containing such particles is administered to the animal or patient prior to the scan. During the data acquisition, the MPI scanner steers the FFP along a deliberately chosen trajectory that traces out the volume of scanning, or at least the field of view. The magnetic nanoparticles within the object experience a changing magnetic field and respond by changing their magnetization. The changing magnetization of the nanoparticles induces a time dependent voltage in each of the receive coils. This voltage is sampled and recorded by a receiver associated with the receive coil. The set of the three sampled and recorded voltages constitutes the acquired data. The parameters that control the details of the data acquisition make up the scan protocol. In the second step of the image generation, referred to as image reconstruction, the image is computed, or reconstructed, from the data acquired in the first step. The image is a discrete 3D amy of data that represents a sampled approximation to the position-dependent concentration of the magnetic nanoparticles in the field of view. The reconstruction is performed by a computer, which executes a suitable computer program. Computer and computer program realize a reconstruction algorithm. The reconstruction algorithm is based on a mathematical model of the data acquisition. As with any reconstructive imaging method, this model can be formulated as an integral operator that acts on the acquired data; the reconstruction algorithm tries to undo, to the extent possible, the action of the model.

Under the reasonable assumptions that the magnetic nanoparticles are confined to the field of view and that the voltages induced by the magnetized nanoparticles depend linearly on the concentration of the particles, the model may be written in the generic form

where t and x denote time and position; V is the field of view; s— (s , s₂. ^)^T is the vector- valued function formed by the recorded signals; G = (G i , G₂, G₃)^T is the vector- valued system function; and c represents the unknown concentration of the magnetic nanoparticles, assumed to be independent of time. (The superscript ^T indicates transposition.) An explicit expression of the system function in terms of the magnetic fields generated by the MPI scanner, the properties of the magnetic nanoparticles, and the properties of the receive coils and receivers has so far been missing. Rather, the system function is determined experimentally in a separate calibration procedure, as described presently. Fourier transformation of (1) with respect to time yields the equivalent equation

where v denotes frequency and s and G(x, ·) are the temporal Fourier transforms of s and G(x, ·), respectively. It has become customary to work with (2) rather than (1). Discretization of (2) leads to a discrete, linear system of equations for the samples of c on a 3D Cartesian grid in the field of view. The reconstruction algorithm computes a regularized, approximate solution of this system of equations. A number of reconstruction algorithms have been devised, and satisfactory images have been obtained. Since the matrix of the linear system of equations is huge and dense, the reconstruction time tends to be long. Working with (1) would not help in this respect.

The data acquisition of MPI is afflicted by a complication: Not only the magnetization of the nanoparticles induces a voltage in the receive coils, but also the applied magnetic field itself. The receive coils pick up the sum of both voltages. For the reconstruction, however, we would like to have the signal generated by the magnetized nanoparticles in isolation. In principle, this wanted signal may be obtained as follows: In a separate data acquisition step, we acquire the signal generated by an empty MPI scanner, which is the signal induced by the applied field. This needs to be done only once for each MPI scanner and scan protocol. Subtracting this separately acquired signal (empty MPI scanner) from the signal at hand (non-empty MPI scanner) leaves the signal induced by the magnetized nanoparticles alone. Some technical effort is required to cope with the enormous dynamic range of these signals. The PCD pursues a different approach and removes the unwanted signal induced by the applied field from the signal at hand by means of a high-pass filter. This is possible, because the unwanted signal is narrowly band-limited. However, the high-pass filter also removes the lower part of the spectrum of the wanted signal and, thus, discards certain information about the concentration of the nanoparticles.

Strictly speaking, the acquired data do not depend linearly on the concentration of the magnetic nanoparticles. Fortunately, the nonlinear effects are very weak. Indeed, good images have been obtained with the generic linear model (1). There seems to be no need to give up linearity.

We turn to the calibration procedure as currently practiced. The purpose of this procedure is to measure the system function. To this end, a point-like probing object with a very small volume VQ and a constant concentration c₀ of magnetic nanoparticles is moved to each grid point x, of the 3D Cartesian grid on which it is desired to reconstruct the image and, while the probing object sits at grid point χ, , the data acquisition is repeated using the scan protocol that has been or will be used for the object proper. The probing object at x, can be modeled by the weighted delta func- tion CQ VQ8(- — Xf ) . Inserting this expression into (1) shows that the calibration step yields the time dependent function G(x_/ , ·), for every grid point x, in the field of view, up to an inessential common scaling factor. This calibration procedure has the following advantages:

• One does not need an explicit mathematical model of the data acquisition of MPI.

• Linearity accepted, a measured system function is not degraded by further simplifying assumptions.

But the current approach also has the following disadvantages:

• The measured system function is degraded by measurement errors. · The resolution of the reconstructed image is limited by the finite size of the probing object.

• Measuring the system function is very time consuming.

• The lack of an explicit representation of the s stem function makes it difficult to reason about the solvability of the reconstruction problem or the achievable resolution.

These disadvantages could be avoided if one had an explicit mathematical model of the data acquisition of MPI that were expressed in terms of the various magnetic fields generated by the MPI scanner, the properties of the magnetic nanoparticles, and the properties of the receive coils and receivers. It might still be necessary to measure some of these quantities, but the experimental effort would be drastically reduced. With luck, an explicit model might also pave the way for faster reconstruction algorithms.

Such an explicit mathematical model is derived herein. This is done in two steps.

The first step is to derive an explicit expression for the system function in the generic model (1) in terms of the various magnetic fields produced by the MPI scanner, the magnetization curve of the magnetic nanoparticles, the sensitivity patterns of the receive coils, and the impulse responses of the receivers; once these ingredients have been measured or determined otherwise, the system function is known, up to an inessential scaling factor. The second step is to relate the explicit generic model derived in the first step with an appropriate integral operator that maps functions defined in the spatial domain and representing particle concentrations onto functions also defined in the spatial domain. We shall see that this operator, M, operates according to the formula

(Mc)(y) = K(x

J f , y)c(x) dx, y e W, (3) v

where W denotes the volume of scanning and the kernel K = {K , K₂, K₃)^T can again be expressed in terms of the various magnetic fields generated by the MPT scanner, the magnetization function of the magnetic nanoparticles, the sensitivity patterns of the receive coils, and the impulse responses of the receivers. In addition, the vector- valued function Mc is related to the vector of the magnetization-induced voltages, v = (vi , v₂, v₃)^T , by v( dt' = (Mc H pit ) ) - ( OWD ). t € [i₈, r_e] , (4) where r_s and r_e are the start and end times of the data acquisition and p : [r_s, t_e]→ W is the FFP trajectory. Thus, integration of the voltage vector v with respect to time yields the function Mc, up to an (unknown) additive constant, sampled along the FFP trajectory p. The voltage vector v can be derived from the signal vector s by deconvolving the impulse responses of the receivers from the respective components of the signal vector. The FFP trajectory can be measured or derived from one's knowledge of the MPI scanner and the scan protocol. Moreover, under favorable yet realistic assumptions, the integral operator M can be turned into a convolution operator. This observation paves the way for fast reconstruction algorithms that make use of 3D uniform and nonuniform fast Fourier transforms. Physically, the kt component of the vector—(Mc)(p(t)) represents the magnetic flux through the A'th receive coil when the FFP is at p(f )■

In some respects, the situation just described is reminiscent of MRI: The data acquisition of 3D MRI has an explicit mathematical model, namely the 3D spatial Fourier transform, . Moreover, the data acquisition of 3D MRI provides a sampled version of the complex-valued function !Fm, where the complex- valued function m represents the transverse magnetization of the object being scanned and the samples are taken along a trajectory in 3D Fourier space that can be derived from the magnetic fields generated by the MRI scanner and from the MRI scan protocol. On the other hand, the action of the operator J⁷ is completely different from that of the operator M, reflecting a radical difference between the imaging principles of MRI and MPI.

2 Modeling the Data Acquisition

In this section, we derive a formal mathematical model of the data acquisition. In doing so, we obtain an explicit representation of the system function that appears in the generic model (1). We also suggest a simplified calibration procedure based on this explicit representation. The underlying MPI scanner is assumed to be the iCD shown Figure 1.

2.1 Preliminaries

Vectors are written as column vectors. The superscript ^T indicates transposition. Points in space will usually be specified by their coordinates in the Cartesian coordinate system shown in Figure 1. A typical point is variously written as x = (x , y, z)^T or as x = (x]_ , x₂, .X3) ^{7 '}, whichever is convenient. The unit vectors along the three coordinate axes are d = (1 , 0, 0)^r, e₂ = (0, 1 , 0)^T , and e₃ = (0, 0, l )^T . The Euclidean norm of a vector a is denoted by || a || . The field of view is represented by a centered subset V C ³ with a simple geometric shape, such as a cube, ball, or cylinder. The volume of scanning is also described by a centered subset W C IR³ with a simple geometric shape. We assume that V C W and will typically have V = W . Data are acquired during the time interval [t^, t_e].

We shall employ the macroscopic version of Maxwell's equations in SI units. For the convenience of the reader, these equations are reproduced here:

V D (5)

3D

V x H (8) In these equations, p and j are the densities of the free charges and currents, and D, E, B, and H are the electric flux density, electric field, magnetic flux density, and magnetic field, respectively. The first of these equations plays no role in the derivation of the model. Inside the MPI scanner, the fourth equation reduces to

V x H = 0. (9)

The relevant equations (6), (7), and (9) are augmented by a constitutive relation between B and H. We assume that this relation takes the form

B = μΗ, (10) where the magnetic permeability μ is a scalar function of position and field strength and describes the magnetic properties of the medium. The magnetic properties of the medium are also described by the magnetization M, which is related to B and H by

B = ο(Η + M). (11)

From (10) and (11),

M = /H (12) where χ = μ/μ ο— 1 is the magnetic susceptibility of the medium. The term "medium" refers to the physical matter inside the MPI scanner during the data acquisition; of this matter, only the magnetic nanoparticles matter in MPI.

Strictly speaking, we need to distinguish between the field H₀ in an empty scanner (without a magnetizable medium in it) and the field H in a non-empty scanner (with a magnetizable medium in it). In an empty MPI scanner, the magnetization is zero so that (6) and (11) imply V · H₀ = 0. Together with (9) we have

V x H₀ = 0, V - H₀ = 0. (13)

In a non-empty MPI scanner, the magnetization is non-zero, so that (6) and ( 11) imply V · H =— V · M. Together with (9) we have

V x H = 0, V - H = -V - M. (14)

We see that H₀ and H are different. An analogous consideration shows that the magnetic flux density B₀ in an empty scanner is also different from the magnetic flux density B in a non-empty scanner. A field like H₀ or B₀ is called an "applied field." A magnetizable medium distorts the applied field. However, in MP1 the distortion is so small that it can be safely neglected. (Including this effect in the model would make it nonlinear.)

A constant current I_\ that flows through an electromagnetic coil generates a position dependent magnetic field Hi. The normalized field H1//1 depends only on the coil and is called the field pattern of this coil. The field pattern can be computed using the Biot-Savart law. A time dependent current / that flows through the same coil generates the magnetic field H1///1. (In MRI, the field pattern is defined as οΗι / I_\ ·) The field pattern of a receive coil is called the sensitivity pattern of this coil.

2.2 The Selection Field

In Cartesian coordinates, the selection field may be written as

H^s(x) = <//*(x). H$(x). llHx))^1'. (15)

The same field can be written in standard cylindrical coordinates as

H^s(r, ,z) = (Η?{Γ,φ,^ζ),Η*(Γ,φ, ^ζ),Η!1(Γ.φ,^ζ))^τ. (16)

Since the selection field is rotationally symmetric about the 2 -axis, its angular component Ηψ is zero, and the two remaining components do not depend on φ. In addition, the Cartesian and cylindrical components of the selection field are related by r =_v.v² + v². (17)

H_x ^s(x) = -H?(r,z). (18) r

H*(x) =¾_r ^s(r,z), (19)

H_z ^s(x) = Il'ir.z). (20)

In an empty MPI scanner, the selection field satisfies the equation V · H^s = 0. Expressing this equation in cylindrical coordinates yields T ^l(^rH^(r,z))) + Ht(r,z) = 0. (21) We now make the idealizing assumption that the z-component H is linear in z and independent of r in some region {(r, φ, z) : 0 < r < r₀, < > £ [0, 2JT), z_niin < z < z_max) . In other words, we stipulate that Hf (r, z ) = Gz . 0 < r < r₀, z_min < z < z_max, (22) for some constant G called the∑- gradient of the selection field. Inserting (22) in (21) shows that the radial component H of the selection field satisfies the family of ordinary differential equations

-I— (rH;^s(r, z)) J =—G, 0 < r < r₀, z_min < z < z_raax. (23)

Solving (23) with the initial conditions H {r, 0) = 0 for \r \ < r₀ gives

ς ¹

H (r, z) = --rG, 0 < r < r₀, z_min < z < z_m;]x. (24) Using (22) and (24) in ( 18)-(20) yields

H^s (x) = G ( - , ζ)^Γ, x² + v² < r₀, z_min < z < z_max. (25)

We have found: The (rotationally symmetric) selection field is linear in x in the region where the z-component is linear in z and independent of r. If (22) holds approximately, (25) will still hold approximately. We assume that (25) holds at least approximately in the volume of scanning, W .

Reasoning differently, equations (13) and the rotational symmetry of the selection field can be shown to imply that the first order approximation H?(r, z) *i -~Hl , z) (26) holds near the z-axis. Taylor expansion of Hf (0, z) around z = 0 yields the first order approximation

H*(r, z) ^ G₀z (27) with G₀ = ·¾#/(0, 0). Using (26) and (27) in (18)-(20) shows that the first order approximation

H^s (x) _¾ Go( - , - , z)^r (28) holds near the isocenter, even when the z-component of the selection field is not linear in z or not independent of r. 2.3 Drive and Focus Fields

The drive field may be written as

3

H^D(x, ∑C_k(x)r£(t), (29) k = l

where /°( is the drive field share of the current flowing through the

pair and Cjfc(x) is the field pattern of this coil pair. Using matrix-vector notation, we rewrite (29) as

H^D(x,f) = C(x)I^D( , (30) where

C(x) := (Ci(x),C₂(x),C₃(x)) (31) is the field matrix and

P(t) = (I?(t),I?(t),lf(t))^T. (32) Similarly, the focus field may be written as

3

H^F(x, =∑C_k{x)I (t), (33) k=i

or more compactly as

H^F(x, = C(x)I^F(i) (34) with

I^F(i) := {I*(t),I*(t),I*(t)y . (35)

For the parallel coil pairs of the iCD, the field pattern Cjt( ) is in V well approximated by (¾e& for some constant c_k. Hence the field matrix C(x) is nearly independent of x G V, diagonally dominant, and certainly invertible.

In the case of the PCD, drive and focus field are generated by different coil pairs. So instead of (30) and (34), we now have

H^D(x, = C^D(x)I^D( , H^F(x, = C^F(x)I^F( , (36) with two different field matrices C^D(x) and C^F(x). The drive and focus fields of the PCD are spatially still modestly homogeneous in V. The matrices C^D(x) and C^F(x) are diagonally dominant and certainly invertible. Remark. The fact that the drive and focus fields of the iCD have identical field matrices is exploited in the derivation of the model in Section 3. The sameness of the field matrices is a consequence of the fact that the iCD, unlike the PCD, has joint coils for the drive and focus fields. If desired, (nearly) identical field matrices within the volume of scanning could also be achieved using disjoint coils.

2.4 The Field Free Point and its Trajectory

The total magnetic field H(x, t ) inside the MPI scanner is the sum of the static selection field and the position and time dependent drive and focus fields,

H(x, t) := H^s(x) + H^D(x, t) + H^F(x, t). (37)

If the drive and focus fields are zero, the total field has precisely one FFP, namely the isocenter. Adding the drive and focus fields shifts the FFP to a different position in the volume of scanning, W. During the data acquisition, drive and focus fields change with time, and the FFP travels along a trajectory p : [i_s, i_c]→ If. This FFP trajectory is implicitly defined by the equation

H(p(i), = 0, t€ [t_s, t_e], (38)

As we shall see in Section 4, the FFP trajectory should be chosen such that it "densely' fills the volume of scanning W, or at least the field of view V C W . The FFP trajectory must also be realizable by the MPI scanner at hand.

It is also possible to fill W with several concatenated or interleaved trajectory segments. In such a case, the data acquisition would be subdivided into several parts, one for each trajectory segment. For notational simplicity, we restrict ourselves to a single segment herein.

The drive and focus fields are added to the selection field by mak ing the currents 1° and Ιζ flow through the appropriate coil pairs. Two problems pose themselves: First, we may want to move the FFP along a desired trajectory p : [f_s, t_e]→ W and need to find the currents I and Ι that do the job. Second, we may be given the currents and need to determine the resulting FFP trajectory. In both cases, we know the selection field and the field patterns of the coils. The first problem is solved by choosing the currents such that (38) is satisfied. To solve the second problem, we need to find, for each t e [i_s, t_t], the point p( that satisfies H(p(f), t) = 0. We work out the details for the iCD. Inserting (30) and (34) into (38) shows that the sum of the currents,

I := I^D + I^F. (39) and the FFP trajectory p are related by

C(p(f ))I( = H^s(p(/ ) ). t e [i₈, i_e] . (40)

If p(r) is given, (40) becomes a linear system of equations for I(t). For each t, this system has a unique solution, which can be easily computed by any linear equation solver. How to split up this solution into I^D( and I^F(f) is left to the designer of the scan protocol. The MPI scanner must be able to realize the chosen currents. When the MPI scanner is driven with the chosen currents, eddy currents in the coils and cross coupling between the coils may cause the actual trajectory to deviate a little from the prescribed trajectory. Conversely, if l(t) is given, (40) constitutes a (generally) nonlinear equation for p(i). This equation has a unique solution, which can be computed by a suitable nonlinear equation solver. Both problems simplify and can be solved analytically, if we assume the ideal selection field H^s(x) = G(— ,— y-, 3)^T and the ideal field matrix C(x) = (e^ , C₂G₂, c₃e₃) .

For later purposes, we introduce a current map I* by the definition

I* (x) := -C(x)-¹H^s(x), x e W. (41)

Physically, I*(x) is the vector of the currents that, when flowing through the drive and focus field coil pairs (in addition to the current dz/ ^s that flows through the z^±- coil), push the FFP from the isocenter to the point x. Mathematically, the I* in (41) is a reparameterization and extension of the I in (40). Indeed, as can be seen from (40) and (41), we have

1(0 = Γ (p( ). (42) The combined drive and focus field can be written as

H^D(x, i) + H^F(x, = C(x)I* (p( ) (43) and the total field becomes

H(x, = H^s(x) + C(x)I*(p( ). (44) With the PCD, (40) needs to be replaced by

C^D(p(/))I^D( + C^F(p(i))I^F( = -H^s(p( ), t e e]- (45)

Given p(t), we can certainly find I^D(0 and l^¥(t) such that (45) is satisfied. However, the sum of these currents is no longer uniquely determined by (45). This makes it possible to design self-intersecting FFP trajectories that satisfy p(? = p(t₂) with ≠ t₂ and I^D(ii)≠ I^O(t₂), l^F( )≠ Ι¾)· In such a case it is impossible to define two current maps I^D* and I^F* such that I^D(i) = I^D*(p(0) and I^F(i) = I^F*(p( )).

To illustrate the above ideas, we design an FFP trajectory for the PCD. Let the data acquisition start at t_H = 0 and end at t_e = T. The first step is to guess a good ideal trajectory for the FFP that also respects the technical capabilities of the PCD. Such an FFP trajectory is the mapping p : [0, T]→ M³ defined by

where the amplitudes ai , a₂ , i3 > 0 are parameters and the frequencies f_\ and f₂ and the duration T are related via some positive integers n _\ , n₂, and «₃ by

_,/i = n₂n₃/ T, f₂ = n iJi ₃/ T. (47)

The mapping p^D(0 ^:= (pi (t), Pi(t)^)^T is a periodic Lissajous curve with period T_L = T / n 3 and frequency ratio

x [—a₂, a₂] x [0, 0] in the plane x₃ = 0. The mapping p^F(0 := (0, 0, p₃(t))^T is a line segment with end points (0, 0,— a₃) and (0, 0, a₃). Note that p(0 = P^D(0 + P^F(0- Our plan is to have the focus field create the line segment p^F(0 while the drive field continuously creates the Lissajous curve p^D(0- Assuming the ideal selection field

H^¾(x) = G(- , -f , x₃y and the ideal field matrices C^u(x) = (c e c%e₂, c$e₃) and C^F(x) = (c ei , c e₂, c^Fe₃) , we find from (45) and (46) that the two fields will serve our plan when the currents are chosen as

I?(t) = ^~a sin(2*/i , /₂ ^D(0 = a₂ sin{2n f₂t), /₃ ^D(0 = 0, (48)

2 c;

G

I[(t) = 0, / (0 = 0, 7^F(0 = ^~a₂(2- - 1). (49) The constants G, can be considered known. The choice of the frequencies f_\ and f₂ in (46) is restricted by the technical capabilities of the PCD, as is the choice of the amplitudes a^. For the PCD, a feasible set of frequencies and amplitudes is f_x = 24639.2 Hz, f₂ = g/i = 25353.6 Hz, a = a₂ = 16 mm, a₃ = 42 mm. The frequency ratio 35/34 makes for a certain density of the FFP trajectory in the xi- and x₂-directions. A commensurate density in ₃-direction is obtained by setting n₃— 84. The Lissajous period is then T_L = n_\ j = ₂f₂ = T/n ₃ = 1 .38048 ms, resulting in a total duration of T = n₃T_L = 0.1 1596 s. During one Lissajous period, the FFP moves by 2a₃/n₃ = 1 mm in ₃-direction. In medical applications of MPI, patient safety imposes a limit of about 40 T/s on the time derivative | ^-B^F(x, ·) | of the focus field B^F := H^F//x₀- In the present case, B^F changes by 2Ga = 0.210 T during T = 0.1 1596 s so that | ^B^F (x. ·) | « 0.210/0.115 T/s, which is well below the safety limit of 40 T/s. The Lissajous curve p^D, the line segment p^F, and the FFP trajectory p = p^D + p^F resulting from these choices of the parameters are illustrated in Figures 5 (a)-(c). The FFP trajectory fills a cuboid of size 32 mm by 32 mm by 84 mm. Preliminary considerations in Section 3 suggest that the density of the trajectory should afford a spatial resolution of the order of 1 mm in this cuboid. We are not yet done, because we would like to fill the whole field of view of the PCD, which is a cube with an edge length of 84 mm. This can be achieved by letting the focus field create a meandering FFP trajectory, as illustrated in Figure 6. The total duration T would increase accordingly, and the drive field would create its Lissajous curve all the time, or at least during the solid (non-dashed) portions of the meandering path in Figure 6.

An FFP trajectory for the iCD could be designed similarly. Of course, the numerical values of the parameters would have to be different.

Currents computed from ideal selection, drive, and focus fields will also be ideal. When such ideal currents are fed into the coils, the resulting real FFP trajectory will deviate a little from the desired ideal trajectory. A modest deviation is acceptable, but it needs to be known for the reconstruction. We can measure the real trajectory or calculate it by solving (38) for p. Alternatively, we may insert the real field patterns, the real selection field, and the ideal trajectory in (40), solve the resulting equation for the currents that would create the ideal trajectory, and feed these currents into the field generating coils of the MPI scanner. The above considerations did not take into account eddy currents in the coils and cross coupling between different coils. These effects could be taken into account by a more sophisticated model.

2.5 Signal Generation and Detection

The applied total magnetic field, call it H₀, orientates the magnetic nanoparticles in the field of view. This results in a magnetization M and a slightly modified total magnetic field H, which together make up the magnetic flux density B = μο(Η + M). The temporal change of the magnetic flux density induces a time dependent voltage Vk in the receive coil pair along the ¾-axis, k = 1, 2, 3. These voltages may be written as

V_k = V_k° + v_k, k = 1 , 2, 3, (50) where V_k° and v_k are the voltages induced by H and by M, respectively. We wish to relate the voltages v_k induced by M with the magnetic nanoparticles causing M.

First, we establish the relation between the magnetic nanoparticles and the magnetization. We define the concentration of the particles at a point x by c(x) = lim (51)

\u\→o \ U \

where U is a 3D region containing x in its interior, \ U \ the volume of U, and n(U) the number of magnetic nanoparticles in U. We assume that c = 0 outside the field of view V . We also assume that the magnetization can be written in the form (cf. (12))

M(x. = (||H(x, || )e(H(x, )c(x), (52) where P : [0, oo]→ R is a magnetization curve with P(0) = 0 and 0 < P'(0) < oo, and e(H) is the unit (or null) vector defined by

The Langevin theory of paramagnetism suggests the magnetization curve

PdH) = mL^_QmH/k_BT), where m is the magnetic moment of the nanoparticles, ΙΪΒ the Boltzmann constant, T the absolute temperature, and L(w) the Langevin function

L(w ) := coth(uj)— l /w (55) illustrated in Figure 7. From the series expansion of coth( ,') one finds

L(w) = -w w H w — ·^■ · (56)

3 45 945

and hence Z (0) = 1/3. For large values of w, L (w) behaves like coth(w), which asymptotically approaches 1. So from (52) and (54), as the strength [|H|| of the magnetic field increases, the strength ||M|| of the magnetization asymptotically approaches a constant, the saturation magnetization. In practice, we will have a mixture of nanoparticles with slightly different magnetic moments, and the Langevin theory will hold only approximately. We can help ourselves by working with something like an "effective" magnetization curve. This curve can be measured and, like the theoretical magnetization curve (54), it will start with a finite slope and asymptotically approach a constant.

Next, we establish the relation between the magnetization and the induced voltages. Fortunately, we can draw on results from MRI. As elaborated in [ 1 , Section 7.2], the magnetization M causes a magnetic flux

<&_k (t) = ₀ / R_fc(x)^rM(x, r) iix, k = 1 , 2, 3, (57)

Jv

through the receive coil along the ¾-axis, where R_fc(x) is the sensitivity pattern of this coil. It suffices to integrate over V in (57), because M = 0 outside V by our assumption that c = 0 outside V . As also shown in [1, Section 7.2], this magnetic flux induces the voltage

Vk it ) = -^ ,ΛΙ) (58) at

across the terminals of this coil. Using matrix- vector notation, we rewrite (57)-(58) as

v( = -μ₀— f R(x)^rM(x, i/x, (59) a t Jv

where

v( = M .

(60) is the vector of the magnetization-induced voltages and

R(x) = (R₁ (x), R₂(x), R₃(x)) (61 ) the sensitivity matrix of the receive coils. With the receive coils of the iCD, the sensitivity matrix R(x) is nearly independent of x, diagonally dominant, and certainly invertible.

Finally, inserting (52) into (59) yields the desired relation between the magnetic nanoparticles and the magnetization-induced voltages, viz

with

K(x, i) = - _oR(x)^r(P(||H(x, f) ||)e(H(x, i))). (63)

We still need to relate the magnetization-induced voltages with the acquired signals. This relationship can be written as

S_k a_k(b_k * V_k) a_k(b_k * (V_k° + v_k)), = 1 , 2, 3, (64) where the time dependent and known function b_k is the impulse response of the receiver associated with the receive coil along the x^-axis, and the unknown constant k is proportional to the gain of this receiver. Assuming properly adjusted gains¹, we may even pretend that the scanner acquires the signals

S_k = a(b_k * (V_k° + ¾)), A: = 1 , 2, 3, (65) where a is a single unknown constant. With some extra effort it is possible to determine , too. We leave it unknown. For the model and the reconstruction, we would like to have the signals s_k = a(bk * Vk), £ = 1 , 2, 3. (66)

' To adjust the gains, put a point-like probing object at the isocenter of the MPI scanner. Probe this object with an FFP trajectory of the form p_k (t) = ^a sin(27rvt)e¾ and measure the signal s_k {t) picked up by the receive coil along the .v* -axis, for A: = 1 , 2, 3. The amplitude a of the FFP trajectory must be so small that the resulting magnetization still depends linearly on the applied field. Adjust the gains of the receivers such that die signals _\ , ¾, and ¾/2 have the same amplitude. We cannot measure these signals directly. However, by repeating the data acquisition with an empty MPI scanner, we can acquire the signals

S^ = a(b_k * V_k°), k = 1 , 2, 3. (67)

(Strictly speaking, the signals thus obtained are induced not by the total field H, but by the applied field H₀. However, as mentioned at the beginning of this section, the difference between H and H₀ is negligible.) This needs to be done only once for each MPI scanner and scan protocol. According to (65)-(67), the wanted signal ,¾ is then found by subtracting S from S_k,

Sk = S_k - Si, k = 1 , 2, 3. (68)

One problem remains: The signals S_k and S have an enormous dynamic range and ordinary A/D converters cannot sample these signals with the required accuracy. To solve this problem, we provide ourselves with a crude approximation S£ to and sample the difference signals S_k— S and — S , which have a lower dynamic range and hence can be sampled using ordinary A D converters. The wanted signal is the difference of the difference signals: s_k = (S_k - ¾) - (S« - Si) , = 1. 2, 3. (69)

Suitable signals S£ can be predetermined by measurement or simulation. Sampled versions of these signals can be stored in the receivers, and analog versions can be generated by D/A conversion of the sampled versions. The acquired data are given by the vector- valued samples s(ti), i = 0, 1 , L, where L > 0 is a large integer and ti = i_s + £ t with Δί = (t_e - t_s)/L, and s( = (s_l (t) , s₂(t) , s₃(t)) ^T . (70)

Unlike the iCD, the PCD acquires the signals (cf. (64))

S™ = a(¾^FCD * V_k) = a(bl^CO * (V_k° + v_k)), k = 1 , 2, 3, (71) where b _c ^CO is a high-pass filter such that b ^CD * V = 0. As a result, the PCD delivers the signals

C^D .PCD

= a(b * v_k), k = 1 . 2, 3. (72) The suppression of V_k° by a high-pass filter is possible, because, with the PCD, V_k° is essentially band-limited to, say, 50 KHz. The useful frequency spectrum of s^^CO ranges up to 1 MHz. However, below 50 KHz, the spectrum of ^^CD is null. This adversely affects the solvability of the reconstruction problem.

2.6 Putting it All Together: The System Function

Combining (60), (62), (66), and (70) yields s(t) = G(x, t)c(x) dx, t e [r_s, i_c] , (73)

with

(ib_x * Gi (x_t .))( \

G(x, 0 = (^ * G₂(x, .))( (74)

\(b₃ * G₃ (x, ·))( /

G(x, -K(x, r) = -^o (x)^7' -( ^>(l|H(x, i) ||)e(H(x, i))). (75)

of dt

A comparison of (73) with ( 1) shows that the function G defined in (74) and (75) is the system function of the iCD, up to the multiplicative constant . The system function of the PCD is obtained by replacing bk in (74) with b ^CD.

2.7 A Simplified Calibration Procedure

From (52), (74), and (75) we see that the system function is determined, up to an inessential multiplicative constant, by the total magnetic field H generated by the MPI scanner, the sensitivity matrix R associated with the receive coils, the magnetization curve P of the magnetic nanoparticles, and the impulse responses bk of the receivers. This observation leads to the following simplified calibration procedure for the determination of the system function, up to a multiplicative constant.

1. Adjust the gains of the receivers so that the received signals satisfy (65). This needs to be done only once.

2. Determine, by simulation or measurement, the sensitivity matrix R(x) for x e V. This needs to be done only once. 3. Determine, by simulation or measurement, the selection field H^s(x) for x £ V . This needs to be done only once.

4. Determine, by simulation or measurement, the drive and focus fields H^D(x, t) and H^F(x, t) for x e V, t e [f_s, i_e]- This needs to be done for each FFP trajectory. At this point, the total field H = H^s + H^D + H^F is known.

5. Measure the effective magnetization curve P of the magnetic nanoparticles to be used. This needs to be done for each new lot of magnetic nanoparticles.

6. Compute the function K(x, i) = - ₀Κ(χ)^Γ (Ρ(||Η(χ, f) ||)e(H(x, t))) for x e V, t _e].

7. Determine, by simulation or measurement, the impulse response bk of the receiver associated with the receive coil along the ^-axis, k = 1 , 2, 3. This needs to be done only once.

8. For each x e V, compute the time derivative G(x, t) := ^ (x, t) and convolve the / -th component of G(x, ·) with the impulse response b - Both operations can conveniently be done in Fourier space.

In the case of the PCD, replace bk with b^FCD in the last two steps of the above procedure.

3 Beautification of the Model

Equations (73)-(75) constitute an explicit mathematical model of the data acquisition of MPI. However, this formulation of the model makes it difficult to answer questions such as how well the acquired signal vector s determines the wanted concentration c and how one can compute this concentration efficiently. In this section, we convert the model (73)-(75) into mathematically more convenient forms that make it easier to answer such questions. 3.1 The Model as an Integral Operator of the First Kind

In this subsection, we seek to relate the model (73)-(75) with an integral equation of the form

g(y) = f K(x. y) f(x) dx, y e W, (76) Jv

where (x) is related with c(x), y with t, K(x, y) with G(x, t), and g(y) with s(t). A formula like (76) gives rise to an associated integral operator, call it X, that maps functions defined on V to vector- valued functions defined on W via the formula

(Xf )(y) :=

J f K(x. y) f(x) dx, y e W. (77) v

Such an integral operator is said to be of the first kind. The adjoint of X is the operator X' defined by

(J 'g)(x) := f K(x, y)^rg(y) < y, x e V, (78) and maps vector- valued functions defined on W to functions defined on V.

The gist of the desired transition from (73) to (76) is how to relate the spatial variable y with the time variable t . In principle, such a relationship is provided by the FFP trajectory, which maps / to p(i). So we would like to reformulate (73) such that the right-hand side depends on t only via p(i). To make this possible, we need to get rid of the convolution with in (74) and of the time derivative in (75).

To get rid of the convolution, we make the idealizing assumption that there exist inverse filters b (t) with b* * b_k = 8, k = 1, 2, 3, (79) where 8 is the Dirac <5-function. This assumption is not overly idealizing: Because the transfer functions of the receivers of the iCD have no zeros in the relevant frequency range, we can find approximate inverse filters b such that b * b^ is a good approximation to the delta function, k = 1 , 2, 3. From (66) and (79) we have b_k* * ¾ = a(b_k* * bk * Vk) = avi_c , k = 1 , 2, 3. (80)

So the deconvolved signal vector v* := (b* * s_{ , b\ * s₂, b * s₃ f (81) v* = «ν. (82)

From (62) we now have v*( ^{= ~} f K( , t)c*(x) dx, t e [r_s, i_c], (83) at J_v with K(x, t) as in (63) and

c* := c. (84)

This last definition also does away with the unknown constant by merging it with the unknown concentration c. The time derivative in (83) can be removed by integration. Indeed, it follows from (83) and the second fundamental theorem of calculus that the deconvolved and integrated signal vector φ(ί) := / ^yy**({tt''))ddtf',, (85) satisfies

Φ(ί) = - / K(x )c*(x)dx+ I K(x,t_s)c*(x)dx, t e [t ,t_e]. (86)

Jv Jv

It remains to make the transition from t to y. At this point, the current map I* introduced in (41) comes in. With its help we define

H*(x, y) := H^s(x) + C(x)I*(y), x e V, y e W. (87)

Physically, H*(x, y) is the total field at x when the FFP is at y. Mathematically, H*(x, y) is an extension and reparameterization of H(x, t). Indeed, comparing (87) with (44) we find that

H(x, = H*(x,p( ), e V, t e [t _e]. (88)

Let us also define the vector-valued function

P(h) := ^o (||h||)e(h), h e ³. (89)

From (63), (88), and (89) we now find that

K(x,r) = -R(x)^rP(H*(x,p( )), XG t€ [i_s,i_e]. (90) So we can rewrite (86) as

Φ (0 =

The right-hand- side of (91) depends on / only via p(t), as wanted. In view of (91), we introduce an integral operator, M, that acts on functions : R³ →· R with support in V according to the formula

(Mf)(y) R(x)^rP(H*(x, y))/(x) </x, (92)

Physically,—(M f )(y) represents the vector of the magnetic fluxes through the three coil pairs caused by the concentration / of nanoparticles when the FFP is at y. Using the operator M, (91) can be rewritten as

Φ (ί) = (Mc* )(p(t)) - («Aic*)(p(/₈)), t [t_s, t_e] . (93)

The second term on the right-hand-side of (93) is an unknown "constant" vector that depends on c* and p(i_s), but not on /. We have found: The deconvolved and integrated signal vector Φ provides the vector- alued function Mc*— (Mc*) p(t_s)), sampled along the FFP trajectory p. If desired, we may rewrite (93) as

Φ( = ^*)(ρ( ), ? e _e]. (94) where M is defined by

(M f y) : = -

This is truly an integral operator of the first kind.

With the PCD, the situation is less satisfactory, for two reasons. First, a transition from t to y via two current maps I^D* and I^F* may not be possible, because such current maps need not exist, as was pointed out in Section 2.4. It may, therefore, not be possible to define an extended and reparameterized total field H*(x, y) by H^s(x) + C^D(x)I^D*(y) + C^F(x)I^F*(y) in analogy with (87). A well defined such field may not even exist, because, in the case of the PCD, the total field at x when the FFP is at y depends on the shifting drive and focus fields individually, and not only on their sum. To help ourselves, we may pretend that C^D(x) = C^F(x) and try to get along with the resulting approximation

H*(x, y) := H^s(x) + C^F(x)I*(y) with ¾ (y) := -C^F(y)^{" 1}H^s(y). (96)

Second, and this is the more serious trouble, the acquired data allow us to compute only something like a high-pass filtered version of Φ . Specifically, if we use the approximation (96), we arrive at

ΦΜ( = (M₀c*)(p(t)) - (M₀c*)(v(t ) - Φ_1Ο (Γ ) , t€ [t_s, (97) where M₀ is obtained from M by replacing H*(x, y) in (92) with H_Q(X, y), and Φΐο( : = Φ— Φ s the unknown complementary low-pass filtered version of Φ .

For mathematical convenience, we extend the matrix-valued functions C and R and the selection field H^s onto ail of ³ in such a way that the extended quantities approximate their ideal counterparts mentioned in Section 2 (i.e., constant matrices, linear selection field). The current map I*, the field H*, and its direction e(H*) are then also defined everywhere on l³. These extensions allow us to integrate in (92) and similar expressions over all of 1R³ and to evaluate M f and similar expressions everywhere in ? . These extensions have no effect on the values of Mf in W , and therefore do not affect the validity of M as a model of the data acquisition of MPT. The precise values of the kernel R(x)P(H*(x, y)) when x V or y ^ W do not matter. We may also multiply the kernel with a window function ω(χ, y) that equals one inside V x W and goes down to zero outside of V x W, without changing the values (M f )(y), y 6 W. Such a modification ensures that the kernel of M and the function M f have a bounded support.

3.2 Related Convolution Operators

An integral operator iven by a formula like

is said to be a convolution operator. The adjoint of X, x e R , (99) is also a convolution operator. We shall show that the operator M defined in (92) is closely related to a convolution operator.

To see this, we insert the expression (41) for the current map I*(y) in (87) to find that

H^*(y, y) = H^s(y) -

= 0. (100) From (100) and (44), the total field satisfies

H* (x, y) = H*(x, y) - H*(y, y) = H^s(x) - H^s(y) + (C(x) - C(y))l^* (y). (101)

The inhomogeneity term (C(x)— C(y))l*(y) comes from the spatial inhomogeneity of the drive and focus fields. This term is small to a high order of ||x— y|| near y and remains small as x departs from y. Neglecting this term in (101) yields the approximation

H* (x, y) ¾ H^s (x) - H^s(y). (102)

Using this approximation in (92) leads to the integral operator

(M₂ f)(y) := / R(x)^rP(H^s(y) - H^s (x))/(x) dx, y e W. (103)

Jv

Since the neglected inhomogeneity term is so small when ||x — y|| is small and P(H*(x, y)) depends only weakly on x and y when ||x — y[| is not so small, the operator M₂ is a highly accurate approximation to the operator M. Replacing M in (93) by M₂ leads to

Φ( = (M₂c^* t)) - (M₂c*)(p( ) +€(t ), t e |/,. /_c | . ( 104) where e(t) is a very small error term that accounts for the difference between M and M₂.

To turn M₂ into something like a convolution operator, we regard the mapping x H H^s(x) as a coordinate transformation and transform the right-hand-side of (103) using the transformation rule for multiple integrals. The result reads

,/)(y) := f R((H^s)^{_ 1} (r))^7'P(H^s(y)— r) i J(H^s>— ¹ (^r) l . ^"((H^s)^{_ 1} (r)) dr, e w, (105) where J _fis-_ri(r) the Jacobian determinant of (i ) ¹, evaluated at r. With the definitions

/(r) := jJ_(Hs)-i(r)|/((H^s) ¹ <r))₅ (106) E(r) ;=R((H^sr¹(r}), (107)

R(r)^rP(s— r) / (r) dr._. s€ H^s { W ) , (108 we have

{M₂f ){i&^>)-¹{$)} ~ ( .₃f)($) s e W(W (1 9) or, equivalent.)^?.

(JW₂/)( )— (¾₃/)(H^s(y)), y ⁽= W (Π0)

Using this identity in (104) elves

Φ{/)— (M._%e*)(p(t)) - (, ₃c*)(p(f_s)) H-€( ),. f e [l„ ^l, wnere

pit) ·= H^s(p(>)), ΓίΓ) := ^((Η*)^-1*»}. (1 i.2) lfH^s is linear as in (25), the coordinate transformation is very simple and

4/G, Alternatively, if H^s is linear or nearly linear, we may replace or approximate

M. by the o erator <A{.₃ given by

The operators <3 and , ₃ defined in (108) and (113) have the same structure, and both look almost like convolution operators, indeed, in the definitions of these operators, we can replace V and W by ³ without affecting their suitabilit as models, We ma even multiply the kernels R(r)^TP(s— r) and R(x)^rP(H^s(y— x» with respective window functions ώίβ - r) and m(y— s) thai equal one in the respective regions

UHW) * H^S( ) := |s - r ; s€ H^s(r), r 6 H^S(K)} (114) and

* F {y— x : y e W, x e V^'} (1 LVS and drop down to zero outside these regions. These modifications ensure that the functions M₃f and M₃f have a bounded support. They are also required for some of the reconstruction algorithms to be presented in Section 4. To highlight the convo- lutional nature of M₃ and M₃, we write down the expression for (M₃f)(s) in terms of its components, including an optional window function OJ(S— r) with ώ(τ) = 1 forr e H^S(W) * H^S(F):

(M₃f)_k(s) = Y [ w(s-T)P_j{s-r)R_jk(T)f(r)dT, sGl³, k = 1,2,3.

(116)

The integrals in (116) are standard 3D convolutions. An analogous expression holds for (M₃f)(y). The adjoint of M₃ is given by

(¾)(r) = f P(s - r)^rR(r)g(s) ds, r e H^S(V). (117) Written out in components and including a window function <S(s— r), this becomes

Analogous expressions hold for (M'₃g)(y). So the adjoints M'₃ and M'₃ are also closely related to convolution operators.

Because the sensitivity matrix appears under the integrals in (108) and (113), the operators M₃ and Mj are not exactly of the form (98). If, however, the sensitivity matrix is independent of x, say R(x) = R₀, then the operators M₃ and M₃ can be rewritten as

(M* f)(s) := Rj / P(s - r) /(r) dr, s e H^S(W), (119)

(M₄f)(y) := Rj f P(H^s(y - x))/(x) dx, y £ W. (120)

Jv

The same modifications that led to (116) turn (R_Q )^~lM₄ and (R_Q )^~1M₄ into proper convolution operators of the form (98).

With the PCD, a reparameterized magnetic field may not exist, and so we may not be able to write down an equation like (101). We can, however, resort to an approximation such as (96) and from then on proceed as above. So the operators M₃, M₃, M₄, and M₄ may perhaps serve as modestly accurate models with the PCD, too. 4 Image Reconstruction

In this section, we state the reconstruction problem, speculate about its solvability, and outline one slow and three fast reconstruction algorithms for solving it. The fast reconstruction algorithms exploit the near convolutional nature of M, the convolution theorem, and the existence of fast algorithms for computing the Fourier transform. Accordingly, we begin this section with an excursion to Fourier space.

4.1 The Fourier Transform and its Relatives

In this subsection, we collect some basic facts about the Fourier transform and its relationship with convolution, the discrete Fourier transform, the fast Fourier transform, and the nonuniform fast Fourier transform. Most of this material is well known. An exception is the nonuniform fast Fourier transform, which is a more recent development.

4.1.1 Fourier Transform

We define the Fourier transform of a complex-valued d -dimensional function / by

= f /(^χ) β^~ί2πχ'ξ ά^χ, i; £ R^d . (121)

We also write for Ψ f ^'. The inverse Fourier transform exists and is given by

(F^→f)(x) = f /(ξ) β^{ί 2π} -^ξ άξ , x £ R^d . (122)

4.1.2 Convolution

The convolution of two d -dimensional functions h and is defined by

( * /My) = f hiy - x)f{x) dx, y€ R^d . (123) The convolution theorem states that h f = h f . (124)

The term deconvolution refers to the problem: Given h and g and knowing that g = h * / , find f . The problem may be tackled in real space and in Fourier space. Solution in real space. Assume that the function k is such that k * h = 8, the Dirac delta function. Then k * g = k *h* f = 8 * / = / . The question is whether such a k exists and how to find it. In principle, by the convolution theorem, k should be such that k = l/h. But since h may have zeros or near zeros and \Η(ξ)\ will certainly tend to zero as \\ξ\\→ oc, trouble is bound to arise: Deconvolution is an ill-posed problem. To overcome the trouble, we need some form of regularization. For example, we may convolve g with a well-behaved function k_K„ that satisfies

A^'rea * /z ¾ 8 or fc_rea;/z ¾ 1. Such a function is k_KZ = I , the inverse Fourier transform of

Λ* = (125)

|/7|² +6²

for some small e > 0. The overbar in (125) denotes complex conjugation.

Solution in Fourier space. Let k := l/h. Then ^^-1 (^ ) = 3^~ (ίΙι ) = 3^^~ι f = ./ · The trouble with the zeros or near zeros of h persists, but the above cure remains applicable too.

4.1.3 Discrete Fourier Transform

The Fourier transform of a function may be approximated by a discrete Fourier transform (DFT):

M

n) = E_mf(x_m) _e-ⁱ²^" , n = 1— ,N. (126) m— l

Likewise, the inverse Fourier transform of / may be approximated by an inverse discrete Fourier transform (iDFT):

N

f(x_m) =∑d_nf( _n)e^{i23t »}-*^>>, m = L.... . (127)

n = l

In these approximate equations,

_ι are two finite sets of grid points in .^d , and {c_m }^₌₁ and {d_n }„_{= l} are sets of suitable weights. If the sampling conditions are right and the weights properly chosen, then the sums in (126) and (127) are good approximations to the integrals in (121) and (122). Specifically, the weights c_m and d„ need to somehow compensate the density of the grid points x_OT in (126) and$_n in (127). The grids {x_m : 1 < m < } and {£„ : 1 < n < N } may or may not be uniform. If both grids are uniform, say with grid spacings Δχ and Δ^, and if M = N and ΑχΑξ = 1 / N, then the two grids are said to be conjugate to each other.

Computing the DFT of iDFT by direct evaluation of (126) or (127) is computationally slow.

4.1.4 Discrete Convolution

The convolution of two functions h and / may be approximated by a discrete convolution:

M

( * f)({y_€) =∑ c_mh(y_t - x„)f(x_m), i = 1 , . . . , L . ( 128)

m = l

In this equation,

is a set of suitable weights. The accuracy of the approximation ( 128) depends on the sampling conditions and on the chosen weights. Direct evaluation of ( 128) can be computationally slow.

4.1.5 Uniform and Nonuniform Fast Fourier Transforms

If the grids {x_m} and {£„} in ( 126) and (127) are uniform and conjugate to each other and if the common grid size M = N satisfies one of several extra conditions, then ( 126) and (127) may be computed efficiently using the Fast Fourier Transform (FFT). A grid size that admits an FFT algorithm will be called FFT-friendly . Actually, there is whole family of FFT algorithms, collectively referred to as "the" FFT. We distinguish between the FFT for the DFT and the iFFT for the iDFT. The FFT/iFFT is a well known classic.

If the grids {x_OT } and {ξ_η} do not support the classical FFT/iFFT, then (126) and (127) may still be computed efficiently using the Nonuniform Fast Fourier Transform (NFFT). Actually, there is a whole family of NFFT algorithms. NFFTs make use of the classic FFT, but are a more recent development. Depending on the type of the input grid and the type of the output grid and on the direction of the DFT, we distinguish between the following types of NFFTs:

• uniform-to-uniform FFT/iFFT (U2U-FFT/iFFT), nonuniform-to-uniform FFT/iFFT (NU2U-FFT/iFFT),

• uniform-to-nonuniform FFT/iFFT (U2NU-FFT/iFFT),

• nonuniform- to-nonuniform FFT/iFFT (NU2NU-FFT/iFFT).

Depending on the grids involved, a U2U -FFT/iFFT may or may not be replaced by a classical FFT/iFFT.

4.1.6 Applications of FFTs and NFFTs

FFTs and NFFTs may be used for the following tasks:

• Discrete Fourier transformation

- Classical case: Uniform input and output grids which are conjugate to each other and have an FFT-friendly size. Best solved by the classical FFT/iFFT.

- Typical nonuniform case: Nonuniform input grid and uniform output grid. Solved by density compensation plus NU2U-FFT/iFFT. This combination is also known as gridding. The term "density compensation" refers to the choice of the weights in (126) or (127). The gridding method is widely used for image reconstruction in MRI [4]. NU2U- and U2NU-FFTs/iFFTs have also been be used for other reconstraction problems [2].

• Resampling

- Typical case: Uniform resampling of a nonuniformly sampled function using gridding followed by a U2U-iFFT or an iFFT.

• Discrete convolution

- Uniform case using the classical FFT: for example, see [5, Section 5.2].

- Nonuniform case using the NFFT: see [3]. 4.2 The Reconstruction Problem

From the acquired data s(¾), i = 0, 1, .... L, we can compute, by discrete de- convolution and integration, the samples (ίι), i = 1, ... , L, of the vector- valued function Φ defined in (85); see Sections 2.5 and 3.1. We also know that there exists an as yet unknown function c* : R³ →· R that is proportional to the position dependent concentration c of the magnetic nanoparticles, has its support in V , and satisfies

Φ(Η) = (Mc*)(v(t_t)) - (Mc*)(j>(t₀)) + « , t = l,...,L, (129) where the ei are unknown, but small error terms. We also know that there exists a well-behaved function g : ³— R³ such that

Φ(ί ) = g(P(¾)), = 1,...^,L. (130) where the p(¾) are the known sampling points of the FFP trajectory p : [f_s, → W . From all this knowledge, we wish to compute a discrete approximation to c* in V. The underlying integral equation

(Mf)(y) - (Mf vi )) = g(y), y e W, (131) is certainly ill-posed [5]. In general, we have (Mf)(y) 0 for y W, i.e. our data are truncated. The following questions arise:

1. Does the integral equation (131) have a unique solution / for each right-hand- side g?

2. How well does an approximation to g determine an approximation to /^'?

3. How can we efficiently compute a discrete approximation to / from a discrete approximation to g?

A sound answer to the first question is not known, but there is reason to believe that the answer is favorable. This optimistic view is motivated by the fact that the operators M₄ and M₄ (with suitable window functions) are invertible. The argument is the same for the two operators; we present it for M₄. We rewrite (M₄ f )(y) as

(M₄f)(y) = Rf f k(y - x)/(x) dx, yell (132)

J L* where

k(z) := ω(ζ)Ρ(Η^δ(ζ)) (133) is the convolution kernel of M₄, modified by a suitable window function ω (cf. Section 3.2). Let g := M₄ f ^'. By the convolution theorem, we have

and hence

k^Tkf = k^Rj T^. (135)

Since k has a bounded support, k is analytic and so has only isolated zeros. Hence is determined everywhere by (135), except at those isolated zeros. But since / has a bounded support too, is also analytic and can therefore be continuously extended at those isolated zeros. In summary, M₄f determines M₄f, which determines , which determines /. We can also write down an explicit, regularized near inverse M\ of M₄. For example, Ml might be defined by

for some small e > 0. We conclude: At least if the term (M /)(p(?_s)) in ( 131 ) were zero, the data g were not truncated, and the sensitivity matrix R(x) = R₀ were constant, we could compute a regularized near solution of (131), by computing a discretized version of M^#g.

Fourier transform theory may also be invoked to reason about the answer to the second of the above three questions. For example, the quicker the Fourier trans- formed kernel k(£) decays as ^→ ∞, the more ill-posed the reconstruction problem is and the less resolution we will get. Also, from the definition of the kernel k and the shape of the magnetization curve P we can conclude that the steeper the slope of P , the slower the decay of k(£ ) as \\ ξ [| → 00. The voxel size of the reconstructed image should reflect the achievable resolution and thus be adapted to the steepness of the slope of the magnetization curve. To be able to reconstruct an image with a desired voxel size, the average and maximum sample distance in the grid formed by the sampling points t(tt) should be commensurate with the voxel size. The third of the above three questions is answered in Section 4.3 below. With the PCD, instead of (129) we have ω(¾) = (Mc* h ) ) - (Mc ^*)(v(to)) - Φ_1ο(ί_£) I = 1 , . . . , L , (137) where Φ|₀(^) is an extra and unknown error term. Strictly speaking, a function with 4>hi(^) = ghi(p(¾')) need not exist when the FFP trajectory intersects itself, but a good approximation to such a function will exist. The temporal high-pass-filtering of the received signals is not exactly equivalent to a spatial high-pass filtering of g. Nevertheless, we can try to solve an equation like

(Mf)(y) - (Mf)(y₀) = g_hi(y), y e W, (138) but the solution will be a sort of strangely high-pass-filtered approximation to the true c*. (The transition from t to p(i) does not exactly turn the temporal high-pass filtering of the measured data s into an equivalent spatial high-pass filtering of the temporally integrated and spatially remapped data g in (130).)

4.3 Reconstruction Algorithms

The standard approach for solving an integral equation like ( 131 ) is to discretize it and to compute a regularized solution of the resulting linear system of equations, using the conjugate gradient method or some other iterative method [5]. Since the matrix of the system is huge and dense, reconstruction algorithms obtained in this way tend to be slow. Variations of this plan arise when we replace the operator M in (131) by one of M^, M₃, M₄, or M₄. The convolutional nature of these latter operators gives the resulting linear systems of equations a special structure, which allows for faster reconstruction algorithms.

4.3.1 Conjugate Gradient Method

Standard discretization techniques applied to (131 ) lead to the linear system of equa- tions

Au - A₀u = d, (139) where u e R is a discrete representation of the unknown function c* on a 3D Cartesian grid in the field of view V;

is the data vector; and

are block matrices whose blocks a! e . 3x _{are cjlosen SLlcn mat a}J_{u js an a}p_ proximation to (Mc^*)(p(ti)). Multiplying the matrix A with a vector is a discrete counterpart of applying the operator M to a function. Similarly, multiplying the transposed matrix A^r with a vector is a discrete counterpart of applying the adjoint M' to a vector- valued function. The matrix A is huge and dense. The sheer size of the system (139) strongly suggests to solve it iteratively. The system (139) is certainly ill-conditioned, reflecting the ill-posed nature of the underlying integral equation (131). In some way or other, regularization [5] has to be built into any solver of ( 139).

A simple iterative method for solving a system like (139) is the Landweber method [5J. The iterative step of this method reads u u w(A - Ao)^r((A - Ao)u - d), (142) where ω > 0 is a relaxation parameter. Thus, the bulk of the work of each iterative step of the Landweber method consists in computing matrix- vector products of the form Ap with p H^M and A^rq with q e IR^3xL. Regularization may be achieved by early termination of the iteration [5]. Even so, the Landweber method is notorious for its slow convergence.

A more attractive approach is to compute the solution of the regularized normal equations

((A - A₀)^r(A - A₀) + cd)u = (A - A₀)^rd, ( 143) where I denotes the identity matrix and > 0 is a regularization parameter. This approach may be seen as an instance of the Tikhonov-Phillips regularization technique [5]. To simplify the notation, we rewrite ( 143) as

Bu = b (144) with

B := (A - Ao (A - Ao) + cd, b := (A— A₀ ) d. (145)

The matrix B is symmetric and positive definite. This makes it possible to solve (144) iteratively with the well proven linear conjugate gradient method [5]. Here is a pseudocode formulation of this method as applied to (144).

Linear Conjugate Gradient Method

Choose an initial guess u and an accuracy threshold€ > 0.

q <— Bu— b

p < - -q

8 - ll Q ll ²

While 8 >€² do:

h ^ B

σ - 8/v^Th

u <- u + σρ

q r- q + ah

^tmp «- II q II ²

8 <- #tmp

The bulk of the work during each iterative step is the computation of Bp, which is in turn dominated by computing matrix-vector products of the form Ap and A^rq. Preconditioning may be used to reduce the number of iterations, but has a cost of its own [5]. It can be shown that the solution of ( 144), and hence the solution of (143), is also the minimizer of the convex quadratic functional

(146) Indeed, the gradient of this functional is Vp(u) = Bu— b. A brief calculation reveals that the functional φ is related to the functional

so φ and p have the same gradient and the same minimizer. Minimizing a functional like (147) is a frequently seen alternative formulation of the Tikhonov-Phillips regularization technique. Generalizations arise when the penalty term a ||u||² in (147) is replaced by other regularizing expressions [5]. Depending on the properties of the penalty term, the resulting generalized functional s may be minimized using the linear conjugate gradient method or some variant of the nonlinear conjugate gradient method [5]. In either case and once again, the bulk of the computational work is computing matrix-vector products of the form Ap and A^rq.

When applied to the problem at hand, the combination of the Tikhonov-Phillips regularization technique and the conjugate gradient method (linear or nonlinear) still suffers from relatively long computation times. To a good deal, this is due to the fact that the matrix A is huge and dense.

4.3.2 Conjugate Gradient Method with NFFTs

Fortunately, by exploiting the near-convolutional nature of the underlying integral operator M, we can dramatically reduce the computational cost of each iterative step of the conjugate gradient method (any variant). To illustrate this, we replace M by M₃ in (129), obtaining

Φ(η) = (M₃c*) ($(n)) - (M₃c*)(p(t₀)) +€_t, = 0, 1 , . . . , L, ( 149)

Since M₃ is a very good approximation to M, the error terms in (149) are slightly different from those in (129), but still small. A suitable window function ώ should be included in the kernel of M₃ (cf. Section 3.2). It is also wise to extend the truncated data from the scanned volume W onto some slightly bigger volume W , letting the extended data go smoothly down to zero in W \ W . Dropping the error terms in (149) and discretizing the resulting equation now leads to a matrix-vector equation of the form

Au - A₀u = d, (150) with d e ^3L, A, A₀ e M^{3l x} , and u e t^ for some L > L and M > M. The convolutional nature of the operators M₃ and M^f ₃ (cf. Section 3.2) means that matrix-vector products like Ap with p R^M and A^rq with q <= R^3L boil down to discrete convolutions (cf. Section 4.2). This fact makes it possible to speed up the evaluation of such matrix-vector products by an order of magnitude using NFFTs, see [3] and the references cited therein. The conjugate gradient method as such need not be changed. If the selection field is sufficiently linear, we may just as well work with M.2 rather than with M₃.

4.3.3 Conjugate Gradient Method with Uniform Resampling and Classical FFTs

As an improvement, we may first resample the nonuniformly sampled data { Φ (t£ ) onto an FFT-friendly uniform grid in W, using gridding followed by an iFFT. From then on, we can use FFTs, rather than NFFTs, during each iterative step. This leads to even faster versions of the conjugate gradient method.

4.3.4 Projections onto Convex Sets

If the sensitivity matrix is nearly independent of x, we can work with M₄ or even with M₄ (if the selection field is sufficiently linear too), of which we know regularized near inverses M^# ₄ and

(cf. Section 4.3). This makes it possible to set up a quickly converging iterative method for solving the associated integral equation. Besides, the computational effort for each iterative step can be kept small owing to the convolutional nature of all these operators. We present the method for M₄, which we rewrite as in (132) and (133); the same game can be played with M₄.

The integral equation to be solved is

(M₄ f )(y) - μ*₄/)(ρ( ) = g(y), y e W. (151)

The suggested procedure is an instance of a general approach known as Projections Onto Convex Sets (POCS) [6]. We ignore discretization issues and present it in continuous terms:

1. Set _,/o = 0 and

g in W,

go (152)

0 outside W.

2. For n = 1 , 2, . . . repeat:

(a) Compute

fn = Mlign-! +

( 153)

Set /„ to zero outside V and clip negative values of „

(b) Stop when „ is good enough.

(c) Compute g in W .

gn (154)

M₄ f„ - (M₄ fn ) (p(i_s) ) outside W .

Being an instance of POCS, this procedure is guaranteed to converge [6] to some limit /*. Moreover, this limit will be "regular" and nearly satisfy (151). As an improvement, one may want to set g₀ not to zero outside W as in (152), but let it smoothly go down to zero outside W . Since the procedure employs a near inverse of M₄ in step 2 (a), it is likely to converge quicker than the conjugate gradient method, which merely employs (a discrete version of) the adjoint of M₄ (represented by A^T or A^T) as a crude approximation to (unless preconditioning is used). Even more, the bulk of the work during the nth iterative step is to evaluate

and M₄ f_n . Both expressions can be computed efficiently in Fourier space, see (134) and (136). As in Section 4.3.3 it is advisable to resample the data first onto a suitable uniform grid; from then on, FFTs can be used to switch between real space and Fourier space. The above method appears to be the most attractive reconstruction algorithm, but it requires that the sensitivity matrix of the receive coils be (nearly) constant. Conclusions

An MPI scanner, nicknamed the Clinical Demonstrator (iCD), has been proposed, featuring

- joint drive and focus field coils with reasonably homogeneous field patterns,

- receive coils with reasonably homogeneous sensitivity patterns,

- receivers with nonzero transfer functions from DC upwards,

- lossless extraction of magnetization-induced signals.

An explicit representation of the system function of the iCD has been derived. This representation involves the various magnetic fields generated by the iCD, the magnetization curve of the magnetic nanoparticles, the sensitivity patterns of the receive coils, and the impulse responses of the receivers.

This explicit representation makes it possible to speed up and simplify the system calibration.

It also leads to an integral operator M such that the deconvolved and integrated acquired signal vector Φ satisfies Φ (ί) = (Mc*) (p(t)) - («Aic*)(p(r_s)) for t e [t , t_c] , where t_s, t_c are the start and end times of the data acquisition, c * is proportional to the concentration of the magnetic nanoparticles, and p :

[4, i_e]→ R³ is the FFP trajectory. An explicit representation of the kernel of M follows from the explicit representation the system function.

The operator M is almost a convolution operator. Several convolution operators that approximate M well have been exhibited.

The near convolutional nature of M makes it possible to discuss the solvability of the reconstruction problem using Fourier transform theory.

The near convolutional nature of M paves the way for fast reconstruction algorithms that make use uniform and nonuniform fast Fourier transforms. Several such algorithms have been indicated. • The joint drive and focus field coils of the iCD give the drive and focus fields the same field patterns. This fact was used to derive the operator M and its relatives. In principle, (nearly) identical field patterns could also be achieved using disjoint coils.

• Unlike the iCD, the Preclinical Demostrator (PCD) under development at Philips Research Europe employs different coils for the drive and focus fields and lacks a lossless signal extraction. Even so, the PCD can benefit from the simplified and speeded-up system calibration. The suggested fast reconstruction algorithms may be less accurate for the PCD, but could still be useful.

References

[1] E. M. Haacke, R. W. Brown, M. R. Thompson, and R. Venkatesan. Magnetic Resonance Imaging: Physical Principles and Sequence Design. Wiley-Liss, New York, 1999.

[2] P. A. Penczek, R. Renka, and H. Schomberg. Gridding-based direct Fourier inversion of the three-dimensional ray transform. JOSA A, 2 1 :499-509, 2004.

[3] D. Potts and G. Steidl. Fast summation at nonequispaced knots by NFFTs. S1AM J. on Sci. Comput., 24:2013-2037, 2003.

[4] H. Schomberg. Notes on direct and gridding-based Fourier inversion methods. In Proc. IEEE International Symposium on Biomedical Imaging: Macro to Nano, pages 645-648, 2002.

[5] C. R. Vogel. Computational, Methods for Inverse Problems. SIAM, Philadelphia, 2002.

[6] D. C. Youla. Mathematical theory of image restoration by the method of convex projections. In H. Stark, editor, Image Recovery: Theory and Application, pages 29-77. Academic Press, 1987. While the invention has been illustrated and described in detail in the drawings and foregoing description, such illustration and description are to be considered illustrative or exemplary and not restrictive; the invention is not limited to the disclosed embodiments. Other variations to the disclosed embodiments can be understood and effected by those skilled in the art in practicing the claimed invention, from a study of the drawings, the disclosure, and the appended claims.

In the claims, the word "comprising" does not exclude other elements or steps, and the indefinite article "a" or "an" does not exclude a plurality. A single element or other unit may fulfill the functions of several items recited in the claims. The mere fact that certain measures are recited in mutually different dependent claims does not indicate that a combination of these measured cannot be used to advantage.

Any reference signs in the claims should not be construed as limiting the scope.

Claims

CLAIMS:

1. An apparatus (100) for forming a concentration image of the concentration of magnetic particles arranged in a field of view (28), which apparatus comprises:

selection means comprising a selection field signal generator unit (110) and selection field elements (116) for generating a magnetic selection field (50) having a pattern in space of its magnetic field strength such that a first sub-zone (52), having a low magnetic field strength and including a field-free point, and a second sub-zone (54) having a higher magnetic field strength are formed in a field of view (28),

drive means comprising a drive field signal generator unit (130) and drive field coils (136a, 136b, 136c) for changing the position in space of the two sub-zones (52, 54) in the field of view (28) by means of a magnetic drive field so that the magnetization of the magnetic particles contained in said object changes locally,

receiving means comprising at least one signal receiving unit (140) and at least one receiving coil (148) for acquiring detection signals, which detection signals depend on the magnetization in the field of view (28), which magnetization is influenced by the change in the position in space of the first and second sub-zone (52, 54), said detection signals being acquired in sampled form while the first sub-zone (52) is moved along a trajectory,

reconstruction means (152) for reconstructing a concentration image of said field of view (28) from said acquired and sampled detection signals (s(ti)) by the steps of: i) determining samples of a vector-valued function (Φ(ίι)) from said sampled detection signal (s(ti)), and

ii) reconstructing the position-dependent concentration (c*) of the magnetic particles in the field of view based on the equation

wherein M is an integral operator modelling the data acquisition and p(ti) are the known positions of the field- free point at the sampling times ti 1 = 0, 1, ..., L.

2. An apparatus (100) as claimed in claim 1,

wherein said reconstruction means (152) is adapted for reconstructing the samples of said vector-valued function (Φ(ί)) from said detection signals (s(ti)) by discrete deconvolution and integration.

3. An apparatus (100) as claimed in claim 1,

wherein said reconstruction means (152) is adapted for reconstructing a concentration image from said samples of the said deconvolved and integrated detection signals by use of a conjugate gradient method.

4. An apparatus (100) as claimed in claim 1,

wherein the integral operator M is defined by

(_.Mf)(y) R(s)^'rP(H*(s_i )) (x) ifx_t

wherein V denotes the field of view, R(x) is the sensitivity matrix of the at least one receiving coil, H*(x, y) is the total magnetic field at position x when the field-free point is at position y, and P(h) is a vector valued function defined by

P(li) := ₀Ρ(|ίι |)«(ί3),

wherein μο is the magnetic permeability of vacuum, P(H) is the magnetization curve and the vector e(h) is the unit vector along the direction of the vector h.

5. An apparatus (100) as claimed in claim 4, wherein the o erator M is replaced by the operator - ^ 3 defined by

(Jk₃/)(s) R(r)^rP(s - r)/(r) <s?r. s e H^s(¥ ),

wherein W denotes the volume of scanning and ' ^{^} R((il ) ί Γ) ) [_s the transformed sensitivity matrix.

6. An apparatus (100) as claimed in claim 4,

wherein the operator M is replaced by the operator Ji₃ defined by (M₃/)( ^') : = j R(x)^rF(H^s(y - x))f{x) dx_t wherein H^s(x) is the magnetic selection field at position x.

7. An apparatus (100) as claimed in claim 4, wherein the operator "'^ is replaced by the operator 4 defined by

(¾/) } := - r) (r) ifr,

wherein Ro is a constant sensitivity matrix.

8. An apparatus (100) as claimed in claim 4, wherein the operator ½₃ is replaced by the operator Jt₄ defined by i_.M₄f)(y) := Rj ^ Y(l (y - x))f(x} dx, wherein Ro is a constant sensitivity matrix.

9. An apparatus (10) as claimed in claim 1,

further comprising focus means (120, 126a-c) for changing the position in space of the said first and second sub-zones of the selection field in the field of view (28) by means of a magnetic focus field.

10. A method for forming a concentration image of the concentration of magnetic particles arranged in a field of view (28), which method comprises the steps of:

generating a magnetic selection field (50) having a pattern in space of its magnetic field strength such that a first sub-zone (52), having a low magnetic field strength and including a field-free point, and a second sub-zone (54) having a higher magnetic field strength are formed in a field of view (28),

- changing the position in space of the two sub-zones (52, 54) in the field of view (28) by means of a magnetic drive field so that the magnetization of the magnetic particles contained in said object changes locally,

acquiring detection signals, which detection signals depend on the magnetization in the field of view (28), which magnetization is influenced by the change in the position in space of the first and second sub-zone (52, 54), said detection signals being acquired in sampled form while the first sub-zone (52) is moved along a trajectory, reconstructing a concentration image of said field of view (28) from said acquired and sampled detection signals (s(ti)) by the steps of:

i) determining samples of a vector-valued function Φ(ίι) from said sampled detection signal (s(ti)), and

ii) reconstructing the position-dependent concentration (c*) of the magnetic particles in the field of view from the equation

wherein M is an integral operator modelling the data acquisition and p(ti) are the known positions of the field- free point at the sampling times t_ls s 1 = 0, 1, ..., L.

11. Computer program comprising program code means for causing a computer to control an arrangement as claimed in claim 1 to carry out the steps of the method as claimed in claim 10 when said computer program is carried out on the computer.