3D tomographic magnetofluorescence imaging of nanodiamonds

: We demonstrate lensless imaging of three-dimensional phantoms of ﬂuorescent nanodiamonds in solution. Magnetoﬂuorescence imaging is employed, which relies on a dependence of the ﬂuorescence yield on the magnetic ﬁeld, and pervading the object with an inhomogeneous magnetic ﬁeld. This ﬁeld provides a ﬁeld-free ﬁeld line, which is rastered through the object. A 3D image of the object is obtained by imaging a set of 2D slices. Each 2D slice image is computed from a set of 1D projections, obtained under diﬀerent projection directions, using a backprojection algorithm. Reconstructed images containing up to 36 × 36 × 8 voxels are obtained. A spatial resolution better than 2 mm is achieved in three dimensions. The approach has the potential for scalability.


Introduction
Imaging the interior of the human body is a tremendously important technique in medical physics. A variety of approaches exist, using X-rays, photons, radio waves, ultrasound, and positrons. Some of these approaches even allow chemical or biofunctional imaging. Detailed chemical imaging is possible in laser microscopy, using fluorescent markers. Such imaging can provide 3D information, with sub-micrometer spatial resolution, but is applicable to small volumes only.
In this paper we pursue the development of a technique [1][2][3] that may in the future lead to biochemical 3D imaging in the animal or human body, in vivo. The method is based on the combination of fluorescent markers distributed within the body, having magnetic-field-dependent fluorescence efficiency, and an inhomogeneous magnetic field pervading the body. The total fluorescence emitted by all markers varies when the field distribution within the body is modified. By implementing an appropriate set of spatial translations of the field as a whole, data can be collected allowing to reconstruct the 3D distribution of markers within the body. This technique represents one variant of "gradient imaging", the imaging class to which also magnetic resonance imaging (MRI) belongs. The method has a potential spatial resolution at the sub-mm level.

Fluorescent nanodiamonds
In this work we employ fluorescent nanodiamonds (FND) as markers. FNDs are nanoscale diamonds containing nitrogen vacancy centers (NVCs), which emit fluorescence when excited with light in the spectral range 460 − 640 nm [4,5]. When FNDs are excited e.g. by 532 nm laser radiation, the emission is in the range 630 − 800 nm [6]. These wavelengths fall into the near-infrared tissue transmission window, in which light has the maximum penetration depth in biological tissue. FNDs appear suitable for in vivo biomedical imaging and for clinical use [5,[7][8][9][10][11][12]. Their high photostability is an advantage compared to typical biomolecules and many common fluorescent markers. The biocompatibility of FNDs and nanodiamonds in general was tested with different animal models and cells, leading to the conclusion that a concentration smaller than 1 mg/ml is generally not toxic [11,12]. The toxicity of higher concentrations depends on the organism model and on the diamond surface purity. FNDs can be produced with diameters as small as 10 nm, which allows them to diffuse into cells, thus enabling their use as markers for The effect used for this imaging technique is the reduction of the microwave-induced fluorescence change (|B|) in absence of a magnetic field. Specifically, ODESR in zero field appears as a fluorescence change by (B = 0) = 0.1% when microwave radiation at appropriate frequency is applied. (|B|) decreases by a factor 2 or more for a field strength |B| ≥ δB = 5 G [3]. This is the effect used in gradient imaging.
Hegyi and Yablonovitch [2,3] demonstrated gradient imaging based on ODESR. Their imaging system combined a permanent magnet QP field, a solenoidal dipole (i.e. spatially homogeneous) field and microwave irradiation at 2.87 GHz. The signal is obtained as the difference of fluorescence when the radiation is on and off. The QP field contained a field-free line with a gradient g 10 G/mm. The combination of this moderate gradient strength and the strong dependence of the fluorescence decrease with magnetic field led to a resolution δρ = 2δB/g 0.8 mm. The field-free line (along the z axis) was scanned by a combination of mechanical and electrical means: by controlling the solenoidal current, the field-free line could be shifted in the y direction, while scanning in x direction was realized by a translation stage moving the object. Imaging also the third spatial dimension (z) was not implemented. The phantoms used were two-dimensional patterns of 2 mm × 2 mm pieces of FND-coated sticky tape, each containing 15 µg of FNDs. For demonstration purposes, these phantoms were also imaged inside a 1 cm × 1 cm × 2 cm piece of chicken breast.

Magnetofluorescence imaging
In the FND ODESR technique described above, the fluorescence change is significant for zero and small field strength. In contrast, the change due to MFE can be significant for field strengths between zero up to a moderate magnitude (few 100 G). Since microwave radiation is absent, and turning the gradient field on and off does not appear to be practical, a modified approach is necessary. Yang and Cohen [1] introduced magnetofluorescence imaging (MFI): here a signal is generated by recording the difference in fluorescence level when the field-free line is slightly displaced in a direction orthogonal to its orientation. This difference signal arises mostly from the fluorescent probes located close to the field-free line.
Yang and Cohen employed a bimolecular system consisting of diamethylanile and pyrene as fluorescent probe. The magnitude of the MFE of this system is dependent on the solvent. This is a disadvantage compared to FNDs, which do not exhibit such an effect, because the NVCs responsible for the fluorescence are mostly located in the bulk. The bimolecular system in the solvent exhibited a decrease of 10% in fluorescence level for field strengths |B| > 90 G. A QP field, generated by four permanent magnets, was employed. Its gradient was substantial, g 260 G/mm, providing a theoretical resolution of δρ = 0.9 mm. Two coils around two opposite magnets allowed spatial shifting and dithering of the field-free line. The samples consisted of 1 mm thick glass pieces of different shapes and did not contain fluorescent particles. For imaging, they were placed into the bimolecular solution inside a container. 2D images were obtained by dithering the field-free line, lock-in detection and scanning the field-free line through the sample with an x-y translation stage. Images of the container with and without the sample were compared to find the samples' boundaries. The technique did not allow to image the samples directly. The authors described the resulting images as "shadow" images.

Modulation MFI
In the present work, we present a 3D tomographic magnetofluorescence imaging setup for FNDs in solution, without the use of microwave irradiation. Our approach, which may be called more specifically "modulation MFI", is a variant of Yang and Cohen's MFI employing modulation of the field-free line. The tomography is implemented by acquiring approximate Radon projections of slices of the sample. A 3D image is obtained from a set of slices. To our knowledge, this non-traditional optical imaging modality has not been shown before, nor has non-traditional 3D imaging of fluorescent markers in macroscopic samples been achieved before.

Magnetofluorescence effect
The key effect at the base of this work, the MFE of NVC, is well known [22,23,25]. To characterize it, we placed a 4 mm dia. × 2 mm cylinder of FND solution on the end face of an electromagnet producing an adjustable magnetic flux B. The solution was irradiated with 100 mW laser light at 532 nm. Fluorescence was collected by a lens and focused on to a solid-state photodetector, yielding a signal S(B). B was measured with a magnetic sensor. Figure 1 shows the fractional fluorescence change MFE(|B|) = S(|B|)/S(B = 0) − 1 vs. |B|. An approximately linear decrease is observed at small to moderate strengths, leveling off beyond 400 G. The half width at half maximum of the response is B c = 200 G. The quantity B c plays a similar role to δB above. We noticed a variation of the details of the MFE function on laser intensity. The fluorescence of FNDs in solution is not expected to show any dependence on the direction of the applied magnetic field.
For imaging of FNDs using a modulation approach or a difference approach, it is a suitable choice to use zero flux strength as the reference. For gradient imaging, the spatial resolution is given by δρ/2 B c /g. In the present work, the gradient at the QP's center has a strength of g 150 G/mm; therefore we expect a spatial resolution of δρ 2.5 mm.

Imaging
We consider a three-dimensional distribution of FNDs described by the number density function n(x, y, z). The imaging task is to obtain this function as faithfully as possible. Figure 2 shows a simplified schematic representation of the measurement procedure. 3D imaging is achieved by combining a set of 2D images of different slices, taken at the positions z 1 , z 2 ,. . . , see Fig. 2(a-e).

2D imaging
2D imaging of the density means obtaining approximately the density in a particular plane, n(x, y, z = z j ). In practice, one seeks to obtain the spatially averaged quantitỹ n(x, y, z j ) = ∫ n(x , y , z )f (x, y, z j ; x , y , z )dx dy dz with the "smearing-out" function f , peaking at (x = x , y = y , z = z ). We may call this a "slice image". Experimentally, one selects by appropriate means the plane z = z j , performs a set of projection measurements therein and then computes a reconstructed quantity that approximatesñ. Indeed, each projection measurement is an approximation of the Radon transform (or projection) of the function n(x, y, z = z j ). The Radon transform [26][27][28][29] is well-known in biomedical imaging, and is used in computer tomography. In the hypothetical ideal configuration, the Radon transform of a 2D slice is a set of line integrals of n in the x-y plane, see Fig. 2(b). In an experimental setting, each line integral is also spatially averaged in the two directions transverse to the projection direction. A set of projections is taken at evenly spaced angles θ i from 0 • to 180 • − ∆θ in steps ∆θ = 180 • /i max . Such a set is shown in a compact representation in Fig. 2(c).
There are different mathematical methods to reconstruct the functionñ(x, y, z j ) (the 2D image) from projections; several available computer codes provide functions for image reconstruction. Here, we use W M : its function InverseRadon uses a filtered backprojection method [30].  yields slice images whose lateral extension does not equal the length of the projection s kmax − s 1 . Instead, their diagonals are approximately equal to the projection's length. The main operation for obtaining a projection measurement consists of orienting the field-free line ì ξ perpendicular to a preselected direction ì s. The orientation of ì s is given by θ i . Different angles θ i are chosen by mechanical rotation of either the object or the magnetic field. During a projection measurement, the field-free line is displaced along ì s, maintaining its orientation θ i , see Fig. 2(b). This displacement is done mechanically in the present work. At each position s k along ì s, one projection data point is taken. Physically, this is a fluorescence level value. The key aspect is that the FNDs located close to the field-free line will contribute to the fluorescence with an amount proportional to their number, and with a slightly higher efficiency due to the magnetofluorescence effect described above. If scattering and absorption effects on both the excitation radiation and on the fluorescence radiation were absent, (thus, in particular, leading to a homogeneous excitation illumination) the increased contribution would be proportional to the line integral of n(x, y, z j ) along the field-free line, averaged over a cylindrical volume around the field-free line with a characteristic, small radial scale.
To obtain a theoretical expression for projection MFI, we approximate the total fluorescence signal emitted into all directions from a point (x, y, z) within the object embedded in a magnetic field B as Here, ζ is the fluorescence efficiency, t 0 is the local illumination level. The total signal detected by a hypothetical detector external to the object and covering the whole solid angle, assuming absence of absorption and scattering of the fluorescence light, is the integral over the whole object, Consider now a QP field with the field-free ì ξ line lying in the plane z j , oriented at an angle θ with respect to the x axis and having the coordinates (x k = s k cos θ + ξ sin θ, y k = −s k sin θ + ξ cos θ).
Here, s k is the smallest distance of the field-free line from the origin x = y = 0, and ξ is the coordinate along it. The QP field strength is given by B also depends parametrically on θ, but this dependence is omitted for clarity. For a fixed θ, we can describe the optical signal using the coordinate system (s, ξ, z) , instead of (x, y, z). The transformation is given by The total signal is given by The integration is over all space, but can be approximated by the line integral of an averaged emitter density, Clearly, the factor MFE describes the spatial averaging of the actual emitter density n. A spatially inhomogeneous illumination t 0 will lead to an image distortion. Only in the ideal case that the illumination t 0 is homogeneous and MFE is a delta function of its argument, would N be proportional to the emitter density n(θ; s k , ξ, z j ) = n(x k , y k , z j ), and the total signal would be a projection along the field-free line ì ξ, This function of s k (θ, z j being parameters) is a Radon transform of the object, the FND density distribution n. Such Radon transforms have to be acquired for a set of angles {θ i } which are input to a backprojection algorithm that outputs an approximation to n(x, y, z j ). The accuracy increases with the number of angles, i max , and of s k values, k max .
In experimental practice, the projection measurement and backprojection will provide the quantity N, Eq. 4, which corresponds to the quantityñ introduced above. It will be affected by additional experimental errors and noise, in addition to the limitations described above.
One important addition to MFI is our variant denoted by modulation MFI of the technique for signal extraction introduced by Yang and Cohen. The position of the field-free line is modulated laterally (in direction of ì s or ì z, perpendicular to the line itself). The ensuing fluorescence modulation amplitude is taken as the signal. Figure 3 shows this principle. The fluorescence modulation arises only from the FNDs close to the field-free line. The fluorescence of the remainder of the object is not modulated (because of the characteristics of the function MFE), and is not recorded in this measurement. There are two modulation options: lateral modulation in the direction ì s or in the direction ì z. In the present work, the modulation amplitudes δz mod and δs mod (Fig. 3) are significantly smaller than the nominal resolution δρ. Then, one data point is approximately given by the derivative where η = z j or s k . In the following, we use the notation X η for ∂ η X. Since neither t 0 nor n depend on η, we can express the integrand as The partial derivative of MFE with respect to η, i.e. z j or s k , can be replaced by the negative of the derivative with respect to z or s, respectively. We denote this variable by η . Integration by parts and taking into account that n(θ; s = ±∞, ξ, z) = 0 leads to This expression is the integrand of S (0) η . Thus, modulation MFI is based on projections of the spatially averaged in-plane (η = s) or out-of-plane (η = z) derivative of the "illuminated FND density" t 0 n. Provided that projections are taken for a sufficiently large number of angles θ i and a dense set of planes {z j } is chosen, the spatial resolution is determined to a large extent by the magnetic gradient g and the derivative of the function MFE, i.e. by δρ. A large signal is obtained when the field-free line moves across an edge of the object in the x-y plane.

3D imaging
A set of slice images taken at different z j are combined to form a 3D image. However, since modulation MFI only captures the derivative of the illuminated FND density, the combination of the slice images yields a 3D image of the illuminated FND density derivative. In the particular case of a homogeneous object (n = const.), z-modulation would result in signals only from the "top" and "bottom" edges of the object, while s-modulation would show the "radial" edges.
Often, one is interested in the FND density itself. For this aim, the projection data can be numerically integrated along the direction of field-free line modulation, and it is on this new data that the projection reconstruction is performed. (The opposite sequence, first reconstruction and then integration, would also work.) Integration renders the homogeneous sectors of an object in the image visible.
The integration is implemented using the trapezoidal rule. In the case of z-modulation, where ∆z = z m+1 − z m is the distance between the planes. In the case of s-modulation an analogous summation is performed.

Simulation
In order to gain insight over the experimental implementation of projection measurements and reconstructions with MFI and how the magnetic field gradient in combination with the MFE influences the resolution, we simulated the imaging of 2D objects, neglecting scattering and absorption, and assuming homogeneous illumination, t 0 (x, y) = t 0 .
Since we are considering a 2D geometry, the magnetic field B(s) (Eq. (1)) can be taken as a linear function. The MFE (Fig. 1) is approximated by a linear function up to B = 450 G. This leads to a triangular function Λ for the MFE with a base width of 2b = 6 mm, see Fig. 4, The scanning of the field-free line along the direction ì s is embodied in this expression. The object is assumed to be a set of circular "patches", centered at the positions (x 0,l , y 0,l ), having radii r l and constant FND densities n 0,l , n l (x, y) = n 0,l circle(x − x 0,l , y − y 0,l , r l ) .
A projection value is computed as  To obtain a complete projection, this expression is evaluated for a set of values {s k }. This is then repeated for a set of angles {θ i }. We choose the step size ∆s = s k+1 − s k , the number of steps k max , and ∆θ to be the same as for the measurements shown in Sec. 5. Fig. 5 presents the simulated imaging of two circular objects of equal density n 0 , having diameters 1 mm and 4 mm. In panel (a) we notice two "bands", one intense and one weak, corresponding to the large and the small object, respectively. In the slice image, panel (b), the reconstructed objects are compared with the position and size of the original objects (circles of solid lines). We find that the maximum values of the reconstructed objects n to be more than a factor 7 larger for the larger object (n 1 = 2.6) than of the smaller object (n 2 = 0.35), even though the densities n 0,l of the objects were the same. The full width at half maximum (FWHM) of the reconstructed smaller object is approximately 3.2 mm and approximately 4.2 mm for the larger one.
To simulate s-modulation, the MFE is expressed as the difference of two triangular functions, displaced by ±δs mod from s k , respectively, The displacement amplitude is chosen as δs mod = 0.175 mm, as in the experiments below. The projections in Fig. 6(a) depict two shifted bands, one of positive and one of negative signal, for each object, instead of one as in the case without modulation. These two bands stem from the edges of the sample. We can apply the reconstruction algorithm to this "raw" data and obtain panel (b). Such a reconstruction does not depict the shape of the sample, since the back transform "interprets" the two shifted bands as two separate objects. By integrating the projections along s, we obtain panel (c). Backprojection yields panel (d). The reconstructed image is very similar to the image in the case of no modulation, Fig. 5(b), which yields a stronger signal than with modulation. This is not applicable to experimental measurements, where obtaining images with modulation yields a stronger and clearer signal than without modulation by employing lock-in detection and filtering of non-modulated signals. However, the reconstructions seen in Fig. 5(b) and Fig. 6(d) show that measurements with modulation should provide the same information after integration as measurements without modulation.

The apparatus
A schematic of the experimental setup is presented in Fig. 7. To create a permanent magnetic field with an approximate quadrupolar character and a well-defined field-free line, four neodymium magnets (DXOXO-N52, KJ M , 2.54 mm dia. × 2.54 mm) were appropriately arranged in an aluminum frame, with the magnets' internal surfaces 40 mm apart. The magnets have a nominal field strength of 6.6 × 10 3 G at their surface. The magnetic field along the s axis was measured using a magnetic field probe, which was placed vertically approximately in-between the two vertical magnets and could be displaced along s (the s axis joins the centers of the two horizontal magnets). Figure 8 shows the obtained magnetic field component B s (s). The gradient in the middle between the magnets was g = 155 G/mm. We compare the measured B s (s) with a theoretical expression B s, theo (s), taken from Refs. [31,32]. With this expression, the QP field along the s axis is approximated as the sum of the magnetic fields of two opposing permanent magnets (see Fig. 8,). On the s axis, the s-component (i.e. horizontal component) of the fields of the two vertical magnets (z direction) can be neglected.
Two coils were placed around two opposing permanent magnets. They were driven with currents in such a way that their fields add constructively in-between the magnets. Modulating the currents then allows to shift the position of the field-free line, along a direction parallel to the line itself. We denote the whole magnet assembly by QPA.  In our setup, we chose to move the total field (QP plus solenoidal), i.e. the QPA, along the z and s directions by mechanical means, while keeping the object fixed. The field-free line therefore always points along the ξ axis. The QPA was mounted under different orientations depending on the type of measurement: the axis of the two solenoids pointed along the z axis when z-modulation was used and along the s axis for s-modulation. Projection angles θ i from 0 • to 180 • − ∆θ are achieved by rotating the sample around the laboratory-fixed axis z in steps ∆θ. Test runs of the apparatus were often performed at a fixed angle θ and for a fixed slice z j , only varying s. This yielded a single projection.
The field-free line was scanned over an s-range of 16 mm with a step size of ∆s = 0.5 mm between measurement points. The coil current was modulated sinusoidally at 240 Hz, with a modulation amplitude of I mod = 5.45 A. This resulted in a measured displacement amplitude δz mod = δs mod = 0.175 mm of the field-free line. The achievable spatial resolution will therefore necessarily be worse than this value. The translation stage displaced the QPA downwards; thus the z 1 = 0 mm plane is located above the FND solution's upper surface (the FND containers are not displaced during the measurement). Samples were placed on a rotary stage within the QPA. They were irradiated from opposite sides with two 532 nm, 100 mW laser diode modules (CW532-100, R L ). Their beams were expanded with engineered diffusers (T ), so that the sample was homogeneously illuminated. The fluorescence light emitted by the sample as well as the excitation light scattered by the sample were focused on a highly sensitive detector (MPPC module C13366-3050GA, H ) with a lens. Two long-pass filters with a cut-on wavelength of 600 nm and one notch filter with a center wavelength of 533 ± 2 nm were placed in front of the detector to suppress the scattered excitation light. The detector signal was sent to a lock-in amplifier (LIA), whose analog output was digitized by a standard data acquisition module. A L program controlled the s-z stage positioning and the sample rotation θ, and recorded the data from the lock-in amplifier for each setting. When not specified otherwise, the LIA was set to an integration time of τ = 300 ms and data was acquired for 1 s per measurement point after a settling time of 1 s before acquisition.
For our experiments, with the particular samples used, the d.c. signal of the detector is typically 2 × 10 3 times larger than the maximum modulation signal.
The LIA signal data is denoted by S in the following. We emphasize that in contrast to the quantity S (0) (and P (0) ) introduced above, S (and P) are affected by absorption and scattering of the fluorescence light. The relationship between S and the FND density is thus more complicated than in the simplified discussion above. In addition, the experimental illumination density is not homogeneous and is not separately determined.

The procedures
The samples consisted of FND solution in small containers, see Fig. 9. The containers were plastic tubes glued to a glass plate. The FND solution used in this work (br100, FND B ) contained 1 mg FND per 1 ml deionized water. According to the manufacturer, the nanodiamonds had a diameter of ∼ 100 nm and were doped with ∼ 900 NVC per single FND. Sample A consisted of two tubes having an inner diameter of 4 mm, liquid fill height ∼ 1.5 mm (object I) and ∼ 1 mm (object II), and placed at a center-to-center distance of ∼ 1.5 mm from each other. Sample B consisted of two tubes of different inner diameter at a distance of 3 mm from each other. One tube had an inner diameter of 4 mm (I), the other of 1 mm (II), both tubes were filled ∼ 2 mm high. The 1 mm tube was a glass tube with an outer diameter of 3 mm.
When placing the samples inside the QPA while preparing for measurements, the objects' long side is aligned with the s axis ( Fig. 9 depicts the objects long side). Hence, 1D projections and the start angle of a set of projections should depict the full diameter of both tubes.

1D imaging
First, 1D measurements of sample A were performed, see Fig. 10. These are projection measurements at a fixed angle θ = θ 1 , that were taken for different z j planes, spaced by ∆z = 0.5 mm. With increasing value z j , the field-free lines shifted from the top to the bottom of the sample.
As shown in Sec. 2.2.1, under simplified assumptions, the recorded signal is approximately given by a weighted mean derivative of the FND density over the area that the field-free line position is modulated. Since the present measurement was performed with z-modulation, only signals from the upper and lower surfaces in z-direction of the sample can be observed. With the particular settings of the LIA, the upper surfaces of the sample produced a positive signal (in the interval z j = 0.5 mm to 3 mm) and the lower surfaces produced a negative signal (interval z j = 4 mm to 6 mm). The projection for z j = 3.5 mm (yellow in Fig. 10(a)) does not show any signal, here the field-free line modulation range lies completely within the objects. Signals can also be observed at the "middle" positions s k = 7 − 7.5 mm, in-between the objects. The fact that they are non-zero might be due to the finite spatial resolution.
To obtain the spatial distribution of the FND solution instead of their derivative, the projections S z were numerically integrated over the modulation direction z to yield P z . The result is shown in Fig. 10(b, c). High values correspond to a higher FND density in the integration (projection) volume perpendicular to the s axis. In panel (a) we recognize a FWHM along the s axis of approximately 3.6 mm, compared with the diameter of the liquid volume of 4 mm. In panels (b, c) the FWHM along the z axis is approximately 2.5 mm for the higher signal and 2 mm for the lower signal. However, the two tubes were filled approximately 1.5 mm and 1 mm high. The discrepancy may be partially due to the employed step size ∆z. The resolution and the FWHM are furthermore influenced by the signal stemming from the area surrounding the field-free line. The FNDs in this area experience low magnetic field and therefore lower magnetic field effects, which contributes to the detected signal, causing a broadening of the signal (as seen in Sec. 3).

2D imaging of dissimilar objects
FNDs are generally considered to be biocompatible in concentrations below 1 mg/ml (see Sec. 1.1). This concentration was used in sample A above. For future imaging in biological tissue with functionalized FNDs, much lower concentrations may occur, determined by the number of functionalized FNDs binding to the targeted tissue. In order to study the concentration dependence of the imaging we performed a 2D projection measurement of a sample C consisting of two 4 mm-diameter tubes, with one tube filled with 20 µl of a 1 mg/ml solution (I), the other with a mixture of 1 µl of the same FND solution and 19 µl deionized water (II). Thus the FND concentration in II was 20 times smaller than in I.
The projections (Fig. 11(a)) show two bands of different strength, as expected. In the resulting reconstruction (panel (b)), the maximum signal of the diluted solution is only smaller by a factor 2.7 compared to the signal of the undiluted solution, while the sum over the pixel values (within black framed areas) differs by a factor of 4.8. These factors are much smaller than the dilution ratio of 20, hence, the relationship between signal height and FND concentration is nonlinear. Considering the variation of the projection signals of the two tubes (panel (a)), we observe that the ratio of the maximum signals changes significantly between the first three angles and the last three angles. This is caused by the positioning of the tubes with respect to the detector. The detector was placed at an angle to the QPA (see Fig. 7), with the objects being closer to the detector for higher values for s k . This could lead to shadowing effects, which are more pronounced when the tube with the undiluted solution blocked the diluted solution from the detectors field of view (see projection angle i = 4).
We explored the spatial resolution in the x-y plane using sample B, measured with z-modulation. Projections were recorded for 20 projection angles from 0 • to 180 • − ∆θ, evenly spaced by ∆θ. Note that fewer angle steps than steps in s direction were used, hence the spatial resolution will not be equal to ∆s. Figure 12 shows the resulting 2D projection and reconstruction. The FWHM of the projections and reconstructed objects are ∼ 3.6 mm and ∼ 2.4 mm for the 4 mm and 1 mm diameter tube respectively. The maximum values of the reconstructed object of the 1 mm tube (II') is ∼ 3.8, significantly higher than for the 4 mm tube (I') ∼ 2.6. The image values summed over area equal to the original object sizes (see Fig. 12(b) framed areas) yield 173 (I') and 31 (II'). One would expect a lower signal from the thinner object (II), as shown in the simulation of Sec. 3, since the absolute number of FNDs within a circle of smaller diameter is lower than within a larger circle. However, the higher signal of the thinner tube was seen over several repetitions of the experiments, so other influences have to be suspected.
When considering the FWHM as an indicator for the spatial resolution, the reconstructions of the 4 mm object would indicate a resolution of 0.5 mm, while the FWHM of the 1 mm object would indicate 1.5 mm. Such a difference in estimated resolution was not expected for different objects that were imaged simultaneously. Since the objects' cross sections are circular, the ξ-integrated FND density is lower on the sides of the objects along s, which results in a lower signal. Hence, the FWHM of a projection through the center of a circular object should be smaller than the actual object diameter, as can be seen in the values for the 4 mm object. This implies that the FWHM is not a good tool to estimate the resolution. However, due to the objects' cylindrical shape, this does not apply for the FWHM of the z dependence. Figure 10(b, c) depicts a FWHM of 2.5 mm and 2 mm in z-direction, while the objects' height was 1.5 mm and 1 mm respectively. This broadening is comparable to the step size ∆z = 0.5 mm and suggests a resolution better than 2 mm. Figure 13(a) depicts 1D projections S s of sample B, obtained using modulation in s-direction at steps of ∆z = 0.5 mm. As described earlier, the projections now emphasize the lateral (in the x-y plane) edges of the sample. In panel (a) this appears as a "dispersive" behavior as function of s. Panel (b) shows the projections integrated along the s direction, P s . We notice a slope in the baselines of these P s . It is in part caused by a small offset in the projections S s (panel (a)) and in part by different absolute values of the maximum and minimum occurring at the edges. A possible origin of the latter effect is the finite step size which leads to an imperfect measurement of S s . Therefore, for the following analysis we correct this effect numerically: for each angle and plane, the sloping baseline is approximated by a line joining the first and last projection data point  and then this value is subtracted from the respective projection. The result of the corrections is shown in panel (c). Figure 14 shows the 2D projection S s for the plane z j = 3 mm performed with s-modulation (panel (a)), the s-integrated projection, P s (panel (b)) and P s after baseline correction (panel (c)). Compared to z-modulation, the s-modulation method shows more artefacts caused by the baseline, which can still be seen after its subtraction (Fig. 14(d)). However, if the given task is imaging a single plane, s-modulation with subsequent integration appears favorable, since it produces a more complete image; specifically, it delineates the objects better, since it is sensitive to edges within the chosen plane. This is not possible with z-modulation. As a result, Fig. 14(d) depicts roughly a cross section of the sample, while Fig. 12 shows the sample's upper surfaces. Another example is the z-modulation projection for z j = 3.5 mm, the yellow line in Fig. 10(a). It is zero, because this slice falls in the middle of the sample. Hence, on its own, this projection does not yield any information. With s-modulation, however, a nonzero signal is obtained, the yellow line in Fig. 13(a), and by integration along s a cut through the 2D slice image is obtained, the yellow line in panel (b).

Comparison of z-and s-modulation mode
The simulations shown in Sec. 3 yield different results for the image values of objects of different diameter compared to the measurements shown here for sand z-modulation, respectively. In simulations with and without modulation the wider object always yields a stronger signal, since a wider object contains a higher total number of FNDs along the ξ axis. We notice that the FWHM of the reconstruction of both the 1 mm and 4 mm object is slightly broader in the simulations than in the experiments; this difference might stem from the simple approximation for the MFE.

3D imaging
We obtained 3D images of sample A by performing 2D z-modulation imaging in a set of planes {z j }. We imaged 8 planes, spaced by ∆z = 1 mm, whose projection data sets S z are shown in Fig. 15. This measurement was performed immediately after the 1D measurement described in Sec. 5.1, therefore the z j planes match and the first angle θ 1 is the same. The projections' main  features are as expected, with two bands (I, II), having different maximum value. Note that the specific band-crossing angle (here i = 7 − 9) depends on the initial orientation θ 1 , while the angles at which the maximum values are attained depend on the location of the detector. The stronger signal of band I originates from the larger object. The projections are not symmetric with respect to the line s k = 8 mm, because the sample was not placed precisely in the middle of the scanned s-interval. We computed the backprojections of the S z for each plane z j and stacked them to form a 3D derivative image, see Fig. 16(a). More artefacts and distortions occur in the slice images, when the signal in a projection measurement occurs near the edges of the reconstructed area. This might cause the elliptical shape of the reconstructed object I' (Fig. 16), whereas the original object had a circular cross section.
The projections S z were integrated over z to obtain P z . As the final step, the backprojection of the P z for each plane z j was computed and stacked to form a 3D image, see Fig. 16(b). Note that the sample appears in only four out of the eight planes measured, and clearly only in two planes (spaced 1 mm apart). The comparison of panels (a) and (b) in Fig. 16 clearly evidences the difference between S z and P z : S z and its backprojection emphasize the borders of the FND distribution along the modulation direction z, whereas the backprojection of P z shows a quantity more directly related to the FND distribution itself. The FWHMs of the projections S z , P z along the s direction and of the slice images, along the x and y direction range from 3.6 mm to 4.5 mm, depending on the z plane. These values are comparable to the diameter of the two objects in sample A, 4 mm.

Discussion and conclusions
Fundamentally, the discussed technique of modulation MFI probes the gradient of the particle density. We have shown that we can probe the gradient both in axial and in radial direction with s-modulation and z-modulation, respectively. This feature can be of interest for some applications. We obtained approximately the particle density by an intermediate step of numerical integration of the data. 3D imaging of different phantoms having different feature size and FND concentration was performed, for the first time. Although the tube samples that were used are not yet a realistic model for future medical applications, our use of FND solution rather than dried FNDs, is already a step in this direction. It was possible to obtain projection sets having 33 × 8 × 14 data points, each covering a "volume" of 0.5 mm × 13 • × 1 mm, and 3D images having 36 × 36 × 8 voxels, each voxel corresponding to a column of approximate volume (0.33 mm) 2 × 1 mm. The required imaging time was 20 min per plane, and 157 min in total. The imaging time of a 2D-measurement with 20 angles, 16 mm scan range with 0.5 mm step size was 28 minutes. A comparison of this imaging duration with previous work was not possible, since no imaging time information was reported [1][2][3].
The FWHM in x and y of projections of the 4 mm-diameter tubes corresponded to the object size, both for sand z-modulation. This suggested a resolution of 0.5 mm. However, the reconstruction of a 1 mm-diameter tube exhibited broadening and suggested a lower resolution of 1 − 1.5 mm. These values should be compared with the naive expectation δρ 2.5 mm. We also found that the imaging of objects with circular cross section complicates understanding of the resolution of the setup. Furthermore, the measurement with diluted and undiluted solution showed that the dependence of the signal on the initial fluorescence intensity is complicated. Clearly absorption and scattering of the fluorescence light, as well as inhomogeneous illumination are important effects. These issues will be studied in future work.

Perspectives
There are several directions into which this work can evolve. A first extension is the implementation of the two modulation modes in an interleaved fashion (after an appropriate extension of the hardware): at every setting (s k , θ i , z j ) of the field-free line, both sand z-modulation are performed in immediate sequence. This increases the total measurement time, but by less than a factor two, since the time overhead for scanning the field-free line across the object is spent only once.
Secondly, the objects imaged here were relatively simple tube samples. Therefore, more complex structures, embedded in scattering media, should be imaged. 3D phantoms could be produced by means of 3D printing.
Thirdly, the artefacts caused by inhomogeneous laser illumination and small detector area should be reduced, which could be done by using distributed light sources and a set of distributed detectors, or by using light guides for illuminating the sample from many directions simultaneously and for gathering the fluorescence light from multiple directions. Since more fluorescence light could then be collected, this would also allow faster imaging or imaging objects with lower FND concentration in similar time.
Better spatial resolution could in principle be achieved with higher magnetic field gradients. Experiments will need to demonstrate this. Further experiments will have to study the influence of scattering media on the imaging performance. Here, an advanced simulation tool for 3D MFI in biological tissue might be helpful for simulating different imaging modalities and experimental configurations before implementing them.
Finally, the most challenging extension is imaging of significantly larger volumes. Whole-body imaging of a 50 cm diameter human body requires a quadrupole magnet with a 1 m bore. This twice larger value is needed for scanning the field-free line across the whole body. If the same gradient strength used here (g 100 G/mm) is to be maintained, the magnetic field at the magnet poles will be of the order of 10 T. Although such values are not impossible, the corresponding magnet will certainly be complex and costly. Possibly, defects other than NVC may exhibit a more favorable magnetic field dependence, thereby reducing the requirements on the magnet. For such large bodies, the MFI approach has potential advantages compared to ODESR imaging. In the latter modality, the larger the volume to image, the higher the power requirement on the microwave source. Furthermore, inhomogeneous distribution of microwave power in the body may occur, leading to image distortions. However, the issue of inhomogeneous irradiation by the laser light is present in both ODESR imaging and MFI imaging and will have to be addressed as mentioned.