1. Introduction

In the last few years fluorescence imaging has become an important tool with many biological and medical applications. Tomographic approaches have been utilized to image novel fluorescent agents with functional and molecular specificity through several millimeters to centimeters of tissue in vivo [1].

Recently developed systems [2, 3] employ non contact detection schemes in order to retrieve the fluorescence measurements. In these systems, the detection of transilluminated light is performed by an optical system that projects the surface of the medium onto a CCD camera producing a large amount of measurement from the pixels of the CCD sensor. Provided that the signal-to-noise ratio is good enough, a larger dataset can provide an improved solution to the inverse problem of finding the fluorescent agent inside a scattering medium. Firstly, the spatial resolution of the reconstruction can be increased using a larger amount of feasible measurements extracted from a large field-of-view [4]. Secondly, the higher information content reduces the illposedness of the problem [5]. Moreover, the use of spectral filters in front of the camera provide multispectral datasets. In the case of multispectral bioluminescence tomography it has been reported that increasing the number of spectral bands renders improved accuracy of the reconstruction results [6, 7, 17]. The use of multispectral reconstructions has also the potential to increase the fluorescent agent contrast by unmixing the fluorochrome signal from the autofluorescence signal [8]. All of the above aspects contribute to very large datasets. In order to solve the inverse tomographic problem the requirements on the computational hardware become immense. In the most common reconstruction schemes the system Jacobian is built and stored in the computer memory. The size of the Jacobian matrix depends on the number of measurements acquired and the resolution of the geometry used for the numerical solution. Practically this imposes a limit for the maximum size of the Jacobian depending on the amount of computer memory that is available. The problem of large datasets has been investigated for the diffuse optical tomography case using an approach based on the analytical solution to the diffusion equation by Wang et. al. [9]. They report on the ability to reconstruct absorption heterogeneities using a number of source-detector pairs in the order of 10⁸.

The approach presented in this report is based on a formulation of the diffusion approximation for the fluorescence case and uses a matrix free formulation. Hence, the explicit calculation and storage of the Jacobian is avoided by replacing it by a vector times matrix operator and a vector times adjoint matrix operator. The method is demonstrated for the reconstruction of fluorescence targets embedded inside scattering medium using simulations and experimental data acquired using a multispectral scheme. We show that the method is beneficial for reconstructions of fluorescence targets using large datasets and multiple wavelengths, since it decreases the computational cost and memory requirements in comparison to the traditional Jacobian methods. Also, the matrix free formulation has the property of requiring equal computer power independent on the number of detectors used, thus making it ideal for imaging detectors.

This paper is organized as follows. In section 2 we describe the problem using a set of two paired diffusion equations for the excitation and the fluorescence wavelength, and construct the derivative operators. Section 3 defines the reconstruction as an optimization problem and section 4 gives an insight to the scalings that were introduced to improve the convergence of the problem. Section 5 presents the matrix-free approach and deals with the implementation issues. We use multiple wavelengths to recover two fluorochromes with different spectral response using numerical simulations in subsection 7.1 and phantom experimental data in 7.2.

2. Formulation of the problem

In multispectral fluorescence reconstruction the fluorochrome concentration distributions of multiple fluorochromes with distinct quantum yield spectra are reconstructed simultaneously from measurements at multiple wavelengths. In analogy to the absorption coefficient [10], we define a wavelength-dependent fluorescence yield coefficient h in a domain Ω:

where η_i(λ) = ε_iγ_i(λ) is the product of the quantum yield γ_i(λ) and the extinction coefficient ε_i for fluorochrome i at wavelength λ, and c_i(r) is the concentration distribution of fluorochrome i.

The forward model for the continuous-wave fluorescence tomography problem is given by the coupled diffusion equations at the excitation and emission wavelengths λ_e and λ_f respectively:

with Robin boundary condition at both wavelengths, [11],

where D and μ _a are the wavelength-dependent diffusion and absorption coefficients, q is a boundary source at the excitation wavelength, ζ is a boundary term incorporating the refractive index mismatch at the surface of the medium, n is the outward surface normal at ξ, and U ^(e) and U ^(f) are the photon density fields at the excitation and emission wavelength, respectively. The contribution of h to the absorption at λ_e is here considered negligible.

The exitance at boundary ∂Ω for both wavelengths is given by the boundary operator

The exitance distribution on the surface is sampled with an array of detectors (e.g. a CCD camera), providing a discrete set of measurements g_d ^(e) and g_d ^(f) at both wavelengths, such that

where w_d is the sensitivity profile of detector d on the boundary.

Assuming that the optical parameters D(r,λ) and μ _a(r,λ) are known, the reconstruction problem consists in finding the fluorochrome concentrations c_i(r). First consider the problem of reconstructing h. Eq. (2–6) define the forward problem which maps distribution h(r,λ_f) to fluorescence data g_d ^(f):

To solve Eq. (2) and (3), we use a finite element model. Parameters D, μ _a and h are expressed as finite-dimensional vectors D, μ _a, h ∈ ℝ^P whose elements are the coefficients of a basis expansion

with basis functions 𝔹 = {b_k(r),k = 1...P}. The fields U are expanded in the basis of a finite element methods (FEM) :

hence the diffusion equation can be written as a linear system

where K ∈ ℝ^N×N is a system matrix assembled from element contributions that depend on the parameter distributions, and U is the vector of basis coefficients representing the photon density distribution in basis expansion 𝕍 = {v_ℓ(r),ℓ = 1...N}. The forward model in the discrete setting consists of solving

with K_λf = K[D(λ3),μ _a(λ_f)], K_λe = K[D(λ_e),μ _a(λ_e)] and ⊙ representing element multiplication. Together with the boundary operator, the forward model can now be expressed in the linear form

where A^(h) is the discrete matrix representation of the forward operator. Since the problem is linear, the Jacobian of the forward operator is the same in matrix representation. Examination of the expression in Eq. (13) indicates that the Jacobian can be constructed using the adjoint formulation

where (sd) denotes a row index constructed from an ordering of source index s and measurement index d, and

are the solutions of the direct and adjoint problems.

3. Inverse problem

In general we pose the inverse problem as an optimisation problem

with regularisation term Ψ and hyperparameter α. Since F(x) = Ax is linear, an iterative solution to Eq. (17) is obtained for example by the damped Gauss-Newton scheme :

where τ denotes the step length. If the regularisation is quadratic, $\Psi \left(\mathbf{x}\right)=\frac{1}{2}{\parallel L\mathbf{x}\parallel }^2$, we have the direct reconstruction formula

where L is for example the Cholesky factorisation of a quadratic Markov Random field. In the results presented in this paper, we consider only the simplest zero-order Tikhonov regularisation L = I.

3.1. Multiple monochromatic reconstruction

The standard way of recovering the fluorochrome concentrations is to take x to be the distributions of h(λ) at the individual fluorescence measurement wavelengths λℓ using

followed by solving

for each basis coefficient k, where it is assumed that the number of wavelengths is greater than or equal to the number of recovered fluorochromes.

3.2. Multispectral reconstruction

In the multispectral problem, instead of reconstructing h, we want to reconstruct the distributions x = {c ₁,c ₂,...,c _N} of multiple fluorochromes c _l simultaneously from data

at multiple wavelengths λ_f. From Eq. (1) we have by chain rule

The complete Jacobian of the problem can then be assembled from blocks for the different wavelengths and chromophores:

where ⊗ denotes the tensor product.

4. Data scalings

Due to uncertainties in laser power, detector gain, and losses we cannot expect the forward model to be compatible with the measured data. To avoid problems due to the model mismatch, we employ data normalization by making use of the available excitation data, and thus we do not require absolute measurements. The optimization problem Eq. (17) is transformed into the rescaled problem

We used two different normalization strategies:

5. Implementation of the matrix-free method

When faced with solving the linear problem Eq. (25), the explicit computation and storage of the matrices A^TA and A^T g _meas ^(f) is costly and often intractable for large scale problems. This is mainly due to the large size of the Jacobian matrix A. As an example, a problem with 30 sources and 475 detectors solved on a geometry with 7812 elements would need about 890MB of storage for each wavelength. The problem is more intense when multiple wavelengths are used, when more samples from the CCD image are required or a high resolution mesh is used to represent the domain in the solver. When using a Krylov solver for the linear problem we require to construct the set of basis vectors

where z = A^T g _meas ^(f) and H = A^TA+αL^TL. In our approach therefore, we represent the forward and adjoint multiplication by the Jacobian implicitly, using a function that returns the result of the operation. To be specific, we use the Matlab ‘gmres’ function to solve Eq. (25), and pass it a function that calculates the forward and adjoint solutions for any intermediate vector generated in the solution iteration. The method is summarised in Algorithm 1.

 

Download Full Size | PDF

6. Materials and methods

We considered a test case where the concentrations of two fluorochromes are to be recovered from measurements at two different wavelengths. We present results from simulated measurements with random Gaussian noise and experimental phantom measurements, showing that the matrix free algorithm is capable of dealing with large data sets reducing the memory and computational costs in respect to the traditional approach.

6.1. Simulation procedure

For the simulation we reproduced the experimental phantom setup of section 7.2 and we assumed a diffusive slab of dimensions 76x69x20 mm ³ containing a rod with concentration 1.7 μM simulating the spectral characteristics of Rhodamine 101 (left) and a rod simulating Rhodamine 6G (right) with fluorochrome concentration 0.4 μM, as shown in Fig. 1(a). In total 30 source positions were used (marked as dots in the back of the figure), and measurements were calculated from 475 positions, placed inside the rectangle area in the front of the figure. Homogenous optical properties were assumed for each wavelength, given in Table 1.

Table 1. Bulk optical properties at the excitation wavelength and emission wavelengths.

View Table | View all tables in this article

 

Fig. 1. (a) A representation of the experimental setup. The dots in the back of the figure denote source positions and the rectangle in the front the image acquired by the CCD camera. The rod on the left contains the Rhodamine 101 while the one on the right the Rhodamine 6G with spectral responses given in the graph on the right. (b) Estimated quantum yield spectra for Rhodamine 101 and Rhodamine 6G.

Download Full Size | PDF

Using the finite element method on a mesh with 6480 voxel elements to solve Eq. (11) we calculated the excitation fields for the wavelength of 532nm and the two fluorescence data set y ⁽⁵⁸⁰⁾ and y ⁽⁶²⁰⁾ for the two wavelengths 580nm and 620nm, respectively. This amounted to 14250 measurements for each wavelength.

Gaussian random noise of 3% was then introduced to the measurements, of both fluorescence and excitation, and the optimisation problem Eq. (19) was solved to recover the concentrations of the two different fluorochromes.

6.2. Experimental setup

Experimental data was acquired using a phantom setup. The bulk scattering media was made from a mixture of water, gelatin, titanium dioxide (TiO ₂) and bovine blood to mimic biological tissue. It was gently stirred to become homogenous and thereafter casted into a slab of size 76x69x20 mm ³. Inside the slab, two cylindrically shaped fluorescent targets with a diameter of 2.5 mm were placed, one containing a 0.4 μM concentration of Rhodamine 6G and the other a 1.7 μM concentration of Rhodamine 101. An estimate of the quantum yield was evaluated by measuring the fluorescence induced spectra in a pure fluorochrome solution. The obtained spectra was normalised with its sum and tabulated quantum yield factor to give the estimated quantum yield spectra. For the two fluorochromes involved, the spectra are presented in Fig. 1(b). The extinction coefficient for Rhodamine 6G and Rhodamine 101 at the excitation wavelength were 23mm ^-1/mM and 6mm ^-1/mM respectively [13]. Non-contact measurements were performed with a CCD-camera (C4742-80-12AG, Hamamatsu) and an objective lens (Nikon f/1.8, focal length 50mm). In front of the lens was a liquid crystal tunable bandpass filter (LCTF VIS 20-35, Varispec) mounted to allow detecting only the fluorescence at λ_f = 580 nm and λ_f = 620 nm, respectively, together with a cut-off filter for blocking the excitation light. A continuous wave (CW) laser (VA-I-N-532, Viasho Laser) emitting at λ_e = 532 nm and generating approximately 10 mW of optical power to the target was used as excitation source. The source was translated along one of the surfaces in a grid pattern by the use of stepper motors. A total of 30 source positions were used together with 1665 detectors sampled across the whole field of view of the CCD-camera, as presented in Fig. 1(a).

The absorption and reduced scattering coefficients were obtained from white-light (Oriel Apex Fiber Illuminator, Newport) transmission measurements with spectrally filtered detection between 540 nm to 720 nm in steps of 20 nm. The spatially resolved transmission data was used to retrieve an effective attenuation coefficient at each spectral band. Time-of-flight spectroscopy was used for assessing the reduced scattering coefficient [14], whereby the absorption coefficient can be computed. The optical properties of interest are presented in Table 1.

 

Fig. 2. Simulated data reconstruction (using 473 detectors per source) Horizontal and Vertical slice along the middle of the slab for the recovered concentrations using simulated data with added 3% noise at the 580nm and 620nm wavelengths for Rhodamine 101 (a) and Rhodamine 6G (b). The computational time for the traditional method of explicit Jacobian was 15minutes 58sec while for the reconstruction using the matrix-free method was 2minutes 44sec.

Download Full Size | PDF

7. Results

7.1. Simulated data

The images of the recovered fluorochromes using the matrix-free method are displayed in Fig. 2 on a horizontal and vertical slice along the middle of the slab. We can see that the location and shape of the target rods were recovered successfully.

For comparison, we reconstructed the same simulated data using the explicit Jacobian method. The resulting images using the matrix-free method are identical to those with the explicit Jacobian. The computational time spend in our 1.8 Ghz machine was about 16 minutes for the traditional method and 2 minutes, 44sec for the matrix-free method.

7.2. Experimental data

The multispectral matrix free algorithm was used and the recovered fluorochrome concentrations are presented in Fig. 3 as horizontal and a vertical slices along the middle of the slab. We notice that the location, the shape and the separation between the two fluorochromes were successfully recovered. The reconstruction took 4 minutes, 32sec. We also tried the same reconstruction using the traditional explicit Jacobian method and the results again were identical to those of Fig. 3. In this case the reconstruction time was was 17 minutes, 5sec.

 

Fig. 3. Experimental data reconstruction using 473 detectors per source Horizontal and vertical slice along the middle of the slab for the recovered concentrations using for phantom experimental measurements at the 580nm and 620nm wavelengths for Rhodamine 101 (a) and Rhodamine 6G (b). This reconstruction using the matrix-free method took 4min. 46sec while the computational time for the traditional method with the explicit Jacobian was 17minutes 4sec.

Download Full Size | PDF

The same reconstruction was performed with the use of a denser sampling on the detectors plane, that resulted to 1665 detectors per source used. The Fig. 4 includes the resulting images from this reconstruction.

There are artifacts appearing close to the detectors plane, in all cases of reconstruction method, especially in the image of the Rhodamine 6G, which we speculate originates from a poor signal-to-noise ratio at the shorter wavelengths.

8. Discussion and conclusions

We have developed and tested a method for the reconstruction of fluorochrome concentrations inside diffusive mediums, taking in respect the latest developments in non contact tomography that allows for a large amount of data to be collected. The proposed method relies on an iterative GMRES solver which calls a functional procedure replacing large matrices with operators that calculate on the run multiplications of matrices with vectors without the necessity for the construction of the matrix. We have shown that the matrix-free method reduces substantially the computational costs and the memory requirements in comparison to the traditional methods that construct an explicit Jacobian matrix. In Table 2 we show the timings and memory requirements for the matrix-free method and the explicit Jacobian for two different measurement setups, one using 1665 × 30 detectors and one with 475 × 30 and for the two wavelengths (580 nm and 620 nm). We should also note that with the matrix-free method the computational cost does not increase when we use more measurements, unlike the traditional methods.

In this paper we have used multispectral data. Increasing the number of spectral bands will inherently create a large number of measurements. The presented method provides a convenient tool for handling such datasets. We believe that this could enhance the fluorochrome contrast for spectrally overlapping fluorphores. This is the focus of our future work.

 

Fig. 4. Experimental data reconstruction using 1665 detectors per source Horizontal and Vertical slice along the middle of the slab for the recovered concentrations using experimental measurements at the 580nm and 620nm wavelengths for Rhodamine 101 (a) and Rhodamine 6G (b). This reconstruction using the matrix-free method took 4min 57sec while the computational time for the traditional method with the explicit Jacobian was 52min 22sec

Download Full Size | PDF

Table 2. Reconstruction times for explicit Jacobian method and the matrix-free using two (580nm and 620nm) wavelengths and two different measurement setups, 475 measurements per source and 1665 positions per source. For the memory allocation calculations, a mesh of 6480 nodes were assumed.

View Table | View all tables in this article