A matrix-free algorithm for multiple wavelength fluorescence tomography

: In the recent years, there has been an increase in applications of non-contact diffusion optical tomography. Especially when the objective is the recovery of ﬂuorescence targets. The non-contact acquisition systems with the use of a CCD-camera produce much denser sampled boundary data sets than ﬁbre-based systems. When model-based reconstruction methods are used, that rely on the inversion of a derivative operator, the large number of measurements poses a challenge since the explicit formulation and storage of the Jacobian matrix could be in general not feasible. This problem is aggravated further in applications, where measurements at multiple wavelengths are used. We present a matrix-free model-based reconstruction method, that addresses the problems of large data sets and reduces the computational cost and memory requirements for the reconstruction. The idea behind the matrix-free method is that information about the Jacobian matrix could be available through matrix times vector products so that the creation and storage of big matrices can be avoided. We tested the method for multiple wavelength ﬂuorescence tomography with simulated and experimental data from phantom experiments, and we found substantial beneﬁts in computational times and memory requirements.


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 8 .
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.

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 Ω: 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: (−∇D(r, λ e )∇ + µ a (r, λ e ))U (e) (r) = q(r) 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 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 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 ∈ R P whose elements are the coefficients of a basis expansion 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 ∈ R 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 = {v ℓ (r), ℓ = 1 . . . N}. The forward model in the discrete setting consists of solving with and ⊙ representing element multiplication. Together with the boundary operator, the forward model can now be expressed in the linear form g 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 J 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.

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, Ψ(x) = 1 2 ||Lx|| 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.

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 λ ℓ usinĝ 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.

Multispectral reconstruction
In the multispectral problem, instead of reconstructing h, we want to reconstruct the distributions x = {c 1 , c 2 , . . . , 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.

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 problemx = arg min We used two different normalization strategies: • The Normalized Born approach [12], which is extensively used in fluorescence tomography.
Where g (e) proj denotes the calculated data in the excitation wavelength. • To get a good balance between the individual fluorescence spectral bands, each spectral band is scaled with its meanḡ ( f ) .

Implementation of the matrix-free method
When faced with solving the linear problem Eq. (25), the explicit computation and storage of the matrices A T A and A T g meas 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 meas and H = A T A + αL T L. In our approach therefore, we represent the forward and adjoint multiplication by the Jacobian implicitly, using a function that returns the result of n · r f n 1 Solve for every source s : end for end for

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.

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 3 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.  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 (580) and y (620) 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.

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 2 ) and bovine blood to mimic biological tissue. It was gently stirred to become homogenous and thereafter casted into a slab of size 76x69x20 mm 3 . 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.

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.

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. . 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.
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.

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 . 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 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.