1. Introduction

Fluorescence Diffuse Optical Tomography (fDOT) aims to visualize and quantitatively characterize molecular events from noninvasive boundary measurements. Fluorescence emissions from specifically designed markers reveal biochemical processes at the molecular level that can be exploited by tomographic imaging. Combining new fluorescent agents with diffuse optical imaging in animal models is likely to become a very powerful tool in the development of new drugs and the study of biochemical processes [1]. In addition, fluorescence lifetime imaging [2–8] is a particularly useful technique because there are many reactions and dynamics of molecular processes that take place on the same time scale as the lifetime of the excited state. The measurement of fluorescence lifetime can provide information concerning the local fluorophore environment in biological tissues. Moreover, protein interactions and conformational changes can be demonstrated using Förster resonant energy transfer, for which fluorescence lifetime provides a reliable read-out [9–13].

In this paper we develop the Discontinuous Galerkin framework (DG) [14, 15] for solving direct and inverse problems in fDOT. We apply this methodology to fluorescence lifetime imaging. DG is the combination of the Finite Volume numerical scheme (FV) and the Finite Elements Method (FEM). The most valuable feature of DG is the ease of mesh adaptation inherited from FV, which is crucial for three-dimensional problems. Secondly, DG method effectively deals with discontinuities in the solution, which may inherently be present in media having, for instance, internal refractive index mismatches [16]. On the other hand, an application of DG numerical scheme to inverse problems has several disadvantages. One of them is the rather high computational cost due to the higher degree of freedom it offers. We show advantages and disadvantages of the DG method by comparing it with a previously developed framework based on FV discretization [17].

The lifetime reconstruction requires time dependent information describing evolution of a physical system. Acquired time dependent data can be Fourier transformed with respect to time to give the equivalent Fourier domain data at multiple harmonic samples. Reconstruction in the Fourier domain has significant advantages over the time domain reconstruction due to its simplicity [7, 17]. Our reconstruction algorithm is designed in the Fourier domain as an iterative solver of a system of equations of Helmholtz type and does not involve full ill-conditioned matrix computations. The algorithm is based on the idea of the reconstruction of a system of parameters iteratively by minimizing the differences of their values estimated for different projection angles and frequencies. The algorithm has been applied to a large experimental dataset acquired by the use of the time gating technique [7, 8, 18–21]. For this type of imaging, the time-gated CCD camera is placed at some distance from the scattering volume, and every pixel of the camera collects photons escaping from the imaging surface within a very short exposure time. Mapping experimental images onto object's surface is not a straight forward procedure. That is because the angular dependence of radiation leaving the surface must be taken into account. In this paper a brief discussion of this problem is also provided.

The paper is organized as follows. Mathematical and computational details of our approach are described in the next section. Section 3 is devoted to experimental details. The last section presents reconstruction results and discussions.

2. Methodology

In contrast to Optical Projection Tomography (OPT), where photons propagate along straight or curved lines, light transport in fDOT is modeled by the diffusion or telegraph equations describing the radiant energy density [22–24]. The diffusive nature of light transport in turbid media such as biological tissues presents significant difficulties in image reconstruction problems and requires development of more sophisticated algorithms than those used in OPT. Thus, the well-known backprojection algorithm cannot be applied in fDOT. Nevertheless, an analogue of the backprojection algorithm for fDOT can be derived by constructing an appropriate cost functional. In this study, the cost functional is built as a sum of excitation and fluorescent error norms and two Lagrangian terms expressing the fact that the energy densities must satisfy the telegraph equation. In addition, the cost functional is penalized by adding one or more regularization terms.

2.1. Inverse problem

We begin by considering a simple experimental setup wherein positions of the laser source and the CCD camera are fixed, but the object under study is rotated over an angle θ. This approach can be easily generalized if necessary. Then, the variational problem is formulated as a minimization problem of the cost functional ℱ:

where the error norm is given as:

The function u_θ and v_θ are the model excitation and fluorescence energy densities, respectively, corresponding to the projection angle θ; the functions e_θ and h_θ are experimentally measured excitation and fluorescence energy densities at the surface of the light scattering object, respectively. The function ξ (θ) is introduced for convenience and for emphasizing similarity with the backprojection operator. It represents the source distribution, which for the case of point sources, as used in this paper, is

where N is the number of source-camera positions. Similarly, the functions χ_θ and ς represent sampling of measurements in space and frequency

where M is the number of discrete points on the imaged phantom's surface; L denotes the number of samples in the Fourier domain (ω); the vector r _θm denotes the surface points visible by the CCD camera corresponding to the excitation source at r _θ. Factors a_m are surface areas around r _θm such that ∫ χ_θ (r)d ³ r gives the total visible area. The form of ℑ_θ is chosen in order to simplify a variational procedure. Thus, the function χ_θ allows to replace a sum over surface points visible by the CCD camera with a volume integral. Analogously, the function ς replaces a sum over samples in the Fourier domain with an integral.

The Lagrangian terms in Eq. (1) are denoted by ℒ_θ and explicitly given by

In expression Eq. (5) 〈·, ·〉 denotes the inner product; the functions ψ_θ and φ_θ are Lagrange multipliers satisfying the same boundary conditions as the energy densities. The Helmholtz operator Λ is given by

where

c is the speed of light in the medium; ρ_θ is the power of the excitation source at r _θ. The unknown reduced scattering and absorption coefficients are denoted as μ′_s and μ_a, respectively. The function q depends on unknown quantum yield, η, the absorption coefficient and lifetime, τ, in the medium as

To avoid a confusion of terminology we term the quantity ημ_a as “fluorescence efficiency”. The quantum yield, η, is understood here as a fraction of two energy densities: the energy density of re-emitted fluorescent photons over the energy density of absorbed excitation ones. The energy loss due to the Stokes shift is already included into this definition.

Note, that unlike the diffusion approximation the telegraph equation assumes dependence of the diffusion coefficient, κ, on the frequency, ω. For low frequencies ∣ω/cμ′_s∣≪1 and μ′_s≫μ_a this dependence is negligible. However, for ∣ω∣ above a few gigahertz and in presence of regions with relatively low values of μ′_s this dependence becomes important.

Next, we define a vector containing unknown optical and fluorescent parameters by x = (μ′_s,μ_a,ημ_a,τ)^T at every point of the domain and choose a dynamic form of the regularization term, ϒ(x), depending on k-th iteration as

where xⁱ_k is i-th component of x _k; and α_i are Tikhonov regularization parameters.

Rather than applying the Gauss-Newton algorithm [25,26] for solving the minimization problem for ℱ we suggest to compute variations of this functional with respect to unknown functions [17, 27, 28]. Setting these variations to zero results in the system of partial differential equations for energy densities, adjoint energy densities (Lagrange multipliers), optical and fluorescence parameters. This system is given below:

where asterisk denotes complex conjugation and functions f_i(θ, x _k) are given by

Note that Eq. (14) contains integration over projection angles θ. Insertion of the function ξ(θ), Eq. (3), into Eq. (14) allows us to rewrite this equation in the form

where xⁱ _0,k = xⁱ_k, and define a subsequence of k-th iteration by letting x _k depend on the projection angle θ_s as

Next, we let the index s run over projection angles in an arbitrary order and associate the index k with samples in the Fourier domain. In this form Eq. (20) presents a variant of the Landweber-Kaczmarz method [29].

Equations (11)–(14) are the Karush-Kuhn-Tucker conditions for optimization of ℑ_θ with conditions ℒ_θ [30]. Equations satisfied by adjoint energy densities, ${\phi }_{{\theta }_n}^{*}$ and ${\psi }_{{\theta }_n}^{*}$ , are Helmholtz equations (Fourier images of corresponding telegraph equations) with source terms placed on the surface of a scattering object visible by the CCD camera. Sources' amplitudes are differences between recorded and computed energy densities. Equations (15)–(18) are used for reconstruction of optical and fluorescent parameters. The physical meaning of these equations can be given by introducing the most probable (averaged) photons' propagation paths. These photons' paths provide the largest contribution to recorded or computed energy densities on the surface of a scattering object from a given source. For instance, for a single source-detector pair averaged photons' propagation paths can be visualized as banana-like distributions connecting source and detector [31, 32]. Thus, parameters are updated iteratively by summing up backprojected differences between recorded and computed energy densities on the visible part of the surface while stepping through projection angles, θ_n, and frequencies, ω_l. Parameters 1/α_i in Eq. (20) are computed at each iteration step as described in [17]. Iterations are terminated when ℑ_θ + ϒ attains its global minimum. Notice, that this approach can be used even for one frequency. For example, for the time independent case (ω = 0) the quantum yield and optical parameters can be reconstructed.

2.2. Implementation

We compute the functions f_i by solving Helmholtz and adjoint Helmholtz equations numerically. The direct solver employs the DG method, which is outlined below. Let us start with the Helmholtz equation written as a system of two first order equations:

Equations (21)–(22) are supplied with boundary conditions according to the approximate intensity

That is, I(s · n<0)=γI(s · n>0) at the open boundary of the scattering domain, where n is the surface normal; s is the unit vector in a particular direction; c q is the energy flux; and γ is a constant depending on the refractive index mismatch [33]. Furthermore, the computational domain is divided into cells, where the solution of Eqs. (21)–(22) is expanded over shape functions ϕ_i(r) satisfying the completeness condition Σ_iϕ_i = 1 inside a cell. For the sake of computational performance we represent optical parameters as piecewise constant functions, following the Finite Volume framework, while all other functions are expanded over piecewise linear basis. Therefore, the energy density and the source term are represented as

Here, and in the rest of this section, a summation over repeated indices is assumed. Multiplying Eq. (21) by ϕ_j and integrating over the cell's volume we arrive at a weak formulation of the direct problem in the local form

where matrix elements are given by

Here V denotes the cell's volume having an outward normal n, and ∂V denotes cell's interface. At this stage the local form, Eq. (25), consists of a system of uncoupled equations. Coupling between cells is provided by replacing q with the interface flux, which is derived below.

As the next step, it is convenient to introduce some useful notations. Let us define jumps and averages across cell's interfaces as

where primes, u′ and q′, denotes corresponding quantities in a neighboring cell. A sum of fluxes at cells' interfaces, Eqs. (26), together with the observation that n = −n′ results in

This expression is further transformed into the following

Noting that the condition [q] · n = 0 is satisfied for the exact solution and taking into account the flux equation q = −κu_i∇ϕ_i Eq. (31) simplifies to

According to this equation we write the term w_j, Eq. (26), in the form

where the interface flux contribution to j-th row and i-th and i′-th columns of the system matrix is

Expression (33) has a simple meaning as a sum of outgoing and incoming fluxes through the shared cell's interface. If the cell interface belongs to the open boundary ∂Ω then boundary conditions are applied

The numerical flux q̂ = −{κu_i∇ϕ_i} in Eq. (32), which replaces q in Eq. (26), is the so-called Bassi-Rebay interface flux [34]. The scheme with the Bassi-Rebay flux is stable for the polynomial interpolation of shape functions ϕ_i of degree higher then 1. For a linear interpolation, the scheme is unstable and must be regularized. For this purpose we impose a set of constraints that the solution is continuous across the cell's interfaces, i.e. [u] = 0, and add the following “zero-term” to the right hand side of Eq. (32):

where parameters β = {−1,0,1} and δ ∈ ℝ⁺ are penalty parameters. In the same way as before we identify v_j in this expression as:

where

These penalty terms improve the condition number of the system matrix.

Some technical details on implementation of DG numerical scheme are outlined below. Firstly, a hexahedral mesh was used in this study. Each computational cell is identified by 8 vertices r _i. Shape functions are chosen in the trilinear form

where {ξ, η, ζ} ∈ [−1, 1]. In order to compute matrix elements we map the cell's interior and interfaces into reference domains, which are a cube and square, respectively. Jacobians of these transformations can be computed by employing quadrature formulae. However, semi-analytical expressions are used in our implementation, which involve computation of the following tensors

plus similar expressions for other interfaces and

where

The octal tree data structure naturally supports mesh adaptation technique. Mesh adaptation by use of this data structure introduces two new cases to consider, namely (i) when a particular cell has many neighbors at the cell's interface, and the case (ii) when the neighboring cell is a coarser cell. The only modification needed to be made to the DG numerical scheme in order to handle these cases is the modification of the incoming flux. That is, the total incoming flux into a cell having many neighbors is just a sum of incoming fluxes from each neighbor. The incoming flux from a coarser neighboring cell is computed as a fraction of an outgoing flux from this neighbor. Technically, it is done by representing a shape function having larger support as a linear combination of shape functions having smaller support on the intersection of their supports.

It is instructive to compare the implementation of the DG method to the Finite Volume (FV) method [17]. FV scheme is similar to the finite difference scheme but can be easily formulated for unstructured meshes. FV scheme can be obtained from DG scheme by replacing u_i and ρ_i in Eq. (25) with their arithmetic averages

and using the completeness condition for shape functions. The numerical flux q becomes 0 inside each cell and, therefore, is defined only on cells' interfaces, where, according to Eq. (22), n · q becomes a δ-function. In the weak form n · q is found in terms of a jump of the energy density ū across the cell's interface [35, 36], see Eq. (29). Both, FV and DG, belong to the family of so-called shock capturing schemes. DG solution space belongs to broken Sobolev's spaces [14] while the solution space of FV is the space of piecewise constant functions. Hence, both methods are well suited for handling discontinuities in parameters as well as in the solution itself [16].

3. Instrumentation

A brief description of the experimental setup is provided in this section. The experimental apparatus consists of three main parts: a pulsed laser, a rotational stage, and a gated camera. The light pulses (pulse width of about 100 ps and central wavelength of 633 nm) are emitted by a diode laser (PDL800, PicoQuant,Germany) operating at the repetition rate of 80 MHz. The light is coupled to a multi-mode graded index fiber and finally focused on the sample placed on a computer controlled rotational stage. The diffusive light, which exits the sample, is imaged by means of an objective (Schneider, Germany) with high numerical aperture (NA = 1.4) on the sensors of the intensified CCD camera (ICCD). In order to separate the fluorescence light from excitation, two filters are used in front of the objective: a long pass glass filter (Schott, Germany), with cutoff wavelength at 665 nm, and a bandpass filter (650–690nm, XF3030 Omega Optical, Battleboro, VT). By removing the filters, the diffusive light at the excitation wavelength is acquired.

 

Fig. 1. (a) Phantom; (b) Surface mesh; (c) Mesh slice showing internal mesh structure.

Download Full Size | PDF

ICCD camera consists of an high repetition rate image intensifier (HRI,Kentech, UK) and a 12-bit Peltier cooled CCD camera (PCO GmbH, Germany). The image intensifier provides a fast gate (approximately 300 ps wide) that allows one to temporally sample the diffused excitation and fluorescence light exiting the phantom. The opening of the gate is synchronized with the master oscillator of the diode laser and delayed by means of a jitter free passive delay generator (minimum temporal step of about 50 ps). Even though a shorter temporal step could be set by the delay generator, this does not provide any advantages due to the jitter of the overall system, which is approximately 30 ps [19]. Since the spatial resolution of the image intensifier is limited to about 128×128 lines, the images are acquired by binning the 1280×1024 CCD sensor to 8×8 pixels. The whole set-up is placed in a light proof black box in order to reduce stray light and secondary reflections that could reach the ICCD sensor at different delays. In order to reduce the acquisition time, the whole measurement procedure is completely automated. For every pixel of the CCD camera the decay profiles of the corresponding temporal signals were fitted by the Green function in infinite space with amplitude, absorption and diffusion coefficients as the 3 fitting parameters. These time-resolved signals were then Fourier-transformed to give complex functions that formed a dataset for reconstruction.

The phantom diagram is shown in Fig. 1.

 

Fig. 2. (a) Image recorded by the CCD camera; (b) corrected image; (c) computed energy density on the surface of a homogeneous cylinder. Intensities of all images are scaled by 10⁵. Images are shown at ω = 0. The source is located on the opposite side of visible part of the surface.

Download Full Size | PDF

It is a solid homogeneous cylinder 40mm in diameter and 100mm in height. The phantom was made of toner and TiO ₂ powder dispersed in the hardener [37]. It has tissue-like optical parameters μ_a = 0.01mm ⁻¹ and μ′_s = 0.83mm ⁻¹, which were measured by a time-resolved spectrophotometer for turbid media based on the Time Correlated Single Photon Counting technique [38]. The phantom has three tubes, two of them were filled with fluorophore (Nile Blue dissolved in methanol) and one was filled with an absorber only. One fluorescent tube is truncated, which makes the inverse problem three-dimensional. All tubes are placed 10mm off the phantom's axis. Tubes A and C contain fluorophore. The concentration of fluorophore in the tube C was four times higher than in the tube A (10⁻⁵ M in the tube C and 2.5 × 10⁻⁶ M in the tube A). However, the quantum yield is the same in both tubes. The tube B was filled with absorber (μ_a = 0.04mm ⁻¹) and its reduced scattering coefficient, μ′_s, has been set to background. Tube A has twice lower value of μ′_s than the background. Tubes A and B have heights 100mm and diameters 4mm. The tube C has smaller diameter 3mm and shorter height 50mm. Its volume is approximately 3.6 times smaller than the volume of the tube A, and, therefore, in spite of higher fluorophore concentration, its brightness is comparable with brightness of the tube A upon excitation. The phantom was probed at three different heights y = {42.5;50.0;67.5} mm. At each height the phantom was rotated by π/6 and imaged. For each camera's position 41 time windows were acquired.

Practically, any reconstruction algorithm requires mapping of experimental images onto the surface of a scattering object. Several mapping procedures were suggested in literature. One of them can be found in works by J. Ripoll et al [39]. Our mapping procedure is different from those and, therefore, we briefly outline it in Appendix. As an illustration of our mapping method, recorded, corrected and computed images for the first projection angle at the excitation wavelength are shown in Fig. 2.

 

Fig. 3. Functions f_i. These functions are used for computing fluorescence and optical parameters. Each row corresponds to a projection angle. Only 3 angles are shown. Each column corresponds to a particular parameter according Eqs. (14)–(18).

Download Full Size | PDF

4. Results and Discussions

As was mentioned above, our reconstruction algorithm, Eqs. (11)–(18), involves repeated solution of four Helmholtz equations. Excitation, fluorescent and corresponding adjoint energy densities are consequently used for computing the functions f_i from Eqs. (15)–(18). An integration of these functions over all projection angles serves as a diffusion analogue of the backprojection operator. However, unlike our previous approach [17], where this backprojection operator was computed for every frequency, we update fluorescent and optical parameters for every projection angle, Eq. (20). This noticeably improves the performance of the algorithm without significant difference in reconstruction results. The drawback of this improvement, however, is that the operator Λ, Eq. (6), must be updated together with optical parameters, i.e. for each projection angle, which was not needed previously. To avoid repeating computation of tensors [Eq. (40)–(43)], we store them in a file. Reconstruction starts with some initial guess where background values of μ′_s and μ_a were chosen. However, any physically meaningful values of the quantum yield and lifetime can serve as the initial guess. Here we set η = 0.001 and τ = 0.01 ns initially everywhere in the computational domain except for the boundary, where η = 0 and τ = 0. Then, all parameters are updated iteratively according to Eq. (20). As an illustration of the algorithm, functions f_i are shown in Fig. 3 for 3 projection angles at ω = 500MHz. Slices are taken at the middle of the phantom.

As we see in Fig. 3, an expected resolution of reconstructed images will not be very high. The lack of resolution is the well-known problem in fDOT in addition to the high computational cost of reconstruction. Computational cost of our algorithm mostly results from repeated solution of Helmholtz equations, while updating parameters does not take much time. The octal tree traversal takes O(nlog₈ n), where n is number of terminal tree nodes, which is much cheaper operation than solving four linear systems.

 

Fig. 4. Reconstruction results. First, second, and third rows show slices at y = 40mm; 50mm; and 60mm respectively. First column shows reconstructed reduced scattering coefficient μ′_s, second column shows the absorption coefficient μ_a, third - the fluorescence efficiency ημ_a, and fourth - the lifetime τ.

Download Full Size | PDF

Reconstruction results are shown in Fig. 4. Each row displays (i) the reduced scattering coefficient μ′_s in mm ⁻¹; (ii) the absorption coefficient μ_a in mm ⁻¹; (iii) the fluorescence efficiency ημ_a in mm ⁻¹, and (iv) the lifetime τ in nanoseconds. Each column displays slices at three different heights y = 40, 50, and 60 mm. Two frequencies were used in reconstruction: 500MHz and 750MHz.

The value of the quantum yield η can be estimated by dividing reconstructed fluorescence efficiency by μ_a. It is seen that the reduced scattering coefficient was reconstructed relatively well. Its minimal value in the tube A is about 0.45mm ⁻¹ at height y = 40mm, which is slightly higher than the true value 0.415mm ⁻¹. The value of μ′_s slightly increases with height achieving 0.52mm ⁻¹ at y = 60mm. As usual, reconstruction artifacts are also present. Reconstruction of the absorption coefficient μ_a is far less accurate. Its value in the tube B is roughly 1.5 lower than it should. Thus, its maximum reconstructed value is 0.027mm ⁻¹ while the true value is 0.04mm ⁻¹. Tubes A and C, which were filled with fluorophore, also appear in reconstruction. This is an expected results due to absorbing properties of fluorophores. Localization of the fluorescent efficiency appears to be relatively good. It is clearly seen that two tubes appear at height y = 40mm and only one at heights 50 and 60mm. The value of the quantum yield is lower than expected. Thus, dividing ημ_a by reconstructed μ_a we obtain approximately [0.18;0.2] at the the center of the tube A, while the true value is about [0.26;0.27]. The lifetime distribution τ, the last column in Fig. 4, has quite high contrast and is well localized. Slices showing the lifetime have background value almost 0. Reconstructed lifetime values are close to the true value of Nile Blue fluorophore, which is 1.2 ns. Thus, (i) at y = 40mm the maximum lifetime is 1.17 ns; (ii) 1.29 ns at y = 50mm; and (iii) 1.14 at y = 60mm. Reconstruction in Fig. 4 shows reasonable errors in optical and fluorescent parameters. The reconstruction error in fDOT must be attributed to the ill-conditioning nature of the inverse problem and can be much higher than it is presented here. It is well-known that it is notoriously difficult to obtain perfectly correct values of parameters in turbid media. The same arguments apply for localization and, especially, for shapes of targets.

 

Fig. 5. Reconstruction results based on the Finite Volume discretization at y = 40mm; 50mm; and 60mm respectively. First column shows reconstructed reduced scattering coefficient μ′_s, second column shows the absorption coefficient μ_a, third - the fluorescence efficiency ημ_a, and fourth - the lifetime τ.

Download Full Size | PDF

4.1. Finite volume method

We compare the DG method to the previously reported FV method [17]. The same octal tree data structure with the same static mesh adaptation technique was employed in reconstruction. FV reconstruction results are shown in Fig. 5.

At first sight, the resolution of reconstructed images is lower than resolution of images computed by DG. It comes as no surprise because FV scheme is only first order accurate in space. Secondly, image mapping adopted in FV is cruder for the same mesh resolution than in DG. That is because DG scheme uses an expansion over shape functions of mapped images while FV computes average values of pixels falling within boundary cells' interfaces. Nevertheless, quantitative accuracy of reconstruction is comparable with DG approach. The reduced scattering coefficient attains values 0.40, 0.32 and 0.42 mm ⁻¹ at the middle of the tube A at height 40, 50 and 60 mm, respectively. The absorption coefficient is almost correct for the tube B at y = 40mm but decreases to 0.03mm ⁻¹ at y = 60mm. Two tubes A and C are well separated at height 40mm and only one tube A is present at heights 50 and 60mm on fluorescence efficiency, as it should. The value of the quantum yield varies in the range [0.20;0.28], which corresponds to reality. Lifetime images does not show separation of tubes A and C at y = 40mm. However, the lifetime value varies in the range [1.0; 1.28] ns, which is a reasonable result.

The main advantage of FV technique is its performance. FV direct solver is at least 8 times faster than DG solver, which obviously speeds up overall reconstruction. That was expected due to different system matrix sizes. Thus, the number of non-zero entries of discretized operator Λ is 64 times larger in the case of DG discretization than in the case of FV. This estimate directly follows from the averaging procedure, Eq. (45). Secondly, the condition number of the operator Λ in FV is smaller, meaning that simple methods for solving linear system can be applied. For instance, due to the diagonal dominance of the operator Λ in the case of FV, methods such as Jacobi or Guass-Seidel iterations can be employed while DG requires more sophisticated algorithms.

Notice that the same optical parameters for excitation and fluorescent light were assumed. Therefore, our reconstructions present somewhat weighted averages of “true” values of optical parameters corresponding to excitation and fluorescent wavelengths. This assumption has the following reason. It is well-known that in absorption-scattering tomography the problem of simultaneous reconstruction of both optical parameters is not unique at a single excitation wavelength. In fluorescence tomography this non-uniqueness is alleviated due to the presence of a fluorescence dataset. Thus, against 4 unknown functions, μ′_s, μ_a, ημ_a and τ, we have 4 distinct datasets: (i) real and imaginary parts of the excitation energy density; (ii) real and imaginary parts of the fluorescent energy density. If we consider distinct optical parameters for excitation and fluorescent light transport, then the non-uniqueness problem will appear again. The same argument applies to the case of multiple lifetimes. Therefore, our assumption is quite reasonable even though values of μ′_s and μ_a can be noticeably different at different wavelengths, as happens in the visible part of the light spectrum. Finally, we would like to emphasize the importance of knowledge of the excitation amplitude, which is a complex number coming from the excitation source term ρ, Eq. (21). As it was suggested previously [17], we use the least squares fit for estimating the absolute value of the amplitude and its phase. Fitting works well when the background optical parameters are known. However, in the case of unknown background values the non-uniqueness problem appears. That is because two additional unknown scalar values are introduced. Thus, starting with any physically meaningful values for background parameters it is always possible to find a reasonably good fit for the excitation amplitude. This obviously affects reconstruction. For instance, higher background value of μ′_s results in increased value of reconstructed ημ_a, while incorrect phase affects reconstruction of μ′_s and τ due to non-linearity of the inverse problem.

In summary, we have demonstrated three-dimensional fluorescence lifetime imaging on the basis of the DG discretization scheme in the Fourier domain with a large experimental dataset. Optical parameters, quantum yield and lifetime were reconstructed simultaneously. The methodology was compared against FV discretization. Our previous variational formulation of the inverse problem was also improved.