A full time-resolved scheme that has been previously applied in diffuse optical tomography is extended to time-domain fluorescence diffuse optical tomography regime, based on a finite-element-finite-time-difference photon diffusion modeling and a Newton-Raphson inversion framework. The merits of using full time-resolved data are twofold: it helps evaluate the intrinsic performance of time-domain mode for improvement of image quality and set up a valuable reference to the assessment of computationally efficient featured-data-based algorithms, and provides a self-normalized implementation to preclude the necessity of the scaling-factor calibration and spectroscopic-feature assessments of the system as well as to overcome the adversity of system instability. We validate the proposed methodology using simulated data, and evaluate its performances of simultaneous recovery of the fluorescent yield and lifetime as well as its superiority to the featured-data one in the fidelity of image reconstruction.
© 2008 Optical Society of America
Near infrared (NIR) fluorescence diffuse optical tomography (fluorescence-DOT) has been extensively explored for contrast enhancement and quantification in optical mammography and brain functional imaging [1, 2]. Combined with the recent advances in red-shifted and NIR molecule-specific fluorescent probes, it is now rapidly evolving as a novel tool to macroscopically visualize and quantify functions of genes and expressions of enzymes and proteins deep inside intact living organisms (mainly small animals), referred to as fluorescence molecular tomography (FMT). This in vivo whole-body imaging modality is a revolutionarily promising and fast growing technique in biomedical research. By binding to or being activated by specific molecules in a biochemical pathway, such as the cell-surface receptor proteins or the enzymes found in tumor cells, the probe detection could be made in the early stage of the target-molecule changes before morphological modification or phenotype becomes apparent, and the related sub-cellular and cellular abnormities be localized and characterized through tomographic imaging of the fluorescence biodistribution. FMT overcomes many disadvantages of the conventional planar photographic modality, such as single projection viewing, strong nonlinear dependency of the signal strength on the depth and the tissue optical properties, as well as low quantification in data interpretation [3–5].
NIR fluorescence-DOT with the continuous-wave (CW) and frequency-domain (FD) modes has been promisingly demonstrated with numerical simulations as well as phantom and in vivo experiments, for both non-specific contrast enhancement and specific molecule tracing [1–5]. The image reconstruction approaches in these two modes are mostly derived from the normalized Born approximation, based on either the relatively simple photon diffusion model [2, 3, 6–10], or the more physically rigorous radiative transfer equation . The CW mode measures only fluorescent intensity as a light of steady intensity passes through tissue. It is of low-cost and achieves a high measuring speed, a high sensitivity and a high spatial resolution (<1 mm) by using a cooled CCD detector for dense boundary sampling [3, 12], but is incapable of measuring lifetime. The FD mode uses an excitation light with its intensity modulated at a radio-frequency, and measures the changes in modulation depth and phase of the fluorescent emission at the same frequency with a heterodyne photomultiplier tube (PMT) system or a gain-modulation image-intensified CCD (ICCD) camera [12, 13]. Although the frequency-domain mode with a single modulation frequency (this is practically the most feasible) can simultaneously reconstruct the fluorescent yield and lifetime images from the intensity and phase signals detected at the boundary, it is unable to analyze multiple fluorescent components, i.e., fluorescent probes of distinct lifetimes targeting at multiple molecular-markers, respectively, and generally requires a modulation frequency of up to 1 GHz and higher for molecular imaging of small animals, which is not technically achievable yet as well as not practical due to exponential decay of the AC amplitude with the frequency [14, 15]. For applications of fluorescence-DOT in molecular imaging, simultaneous recovery of the fluorescent yield and lifetime distributions provides information concerning not only the location but also the local environment of the fluorophore (oxygen, [Ca+] and pH etc.), while analysis of multiple molecular-markers allows assessment of multi-gene controlling mechanism in a disease progression. To attain all these worthy performances without resorting to technically complicated spectrum-resolved methods, it is preferable from both the theoretical and instrumental points of view that the imaging modality be implemented in a time-domain mode [7–9], which employs an ultra-short (femtosecond to picosecond) pulsed light for the excitation and detects the time-dependent fluorescent light using a time-resolved instrument, such as a time-gated ICCD camera [14, 15] or a time-correlated single photon counting (TCSPC) system . By the inverse Fourier Transform, the time-domain system is mathematically equivalent to a frequency-domain one employing a wide range of modulation frequency and therefore possesses the richest information about the optical and fluorescent properties of the chromophores and fluorophores inside a tissue volume. In addition, making efficient use of the time-resolved data would significantly improve fidelity of the image reconstruction, and therefore truly disclose the interactions and processes in complex biological systems.
Analogous to the time-domain DOT, the time-domain fluorescence-DOT can also be performed in either a featured-data scheme [8, 9] or a full time-resolved one, with the former presenting computational simplicity and redundancy reduction while the latter providing better quantification of the reconstruction but at higher computation cost [7, 17]. Nevertheless, the full time-resolved scheme is technically favorable because it helps set-up a reliable reference to the evaluation of the other featured-data or intensity-independent modes . More importantly, it is practically more feasible in molecular imaging than in ordinary DOT since a self-normalized scheme might be introduced not only to avoid the required scaling-factor calibration and the spectroscopic characterization of the system that are usually onerous and difficult to accurately achieve in the Born normalization of CW- and FD-fluorescence-DOT, but also to effectively reduce the adversity of system fluctuations (both from the source power and detector instability) that may occurs between excitation and emission measurements. In principle, the full time-resolved scheme employs an efficient forward model to predict the time-resolved fluorescent flux at the boundary from an initial set of the fluorescence parameters and then updates them to achieve an optimal fitting between the predicted and measured data-sets. In this work, a hybrid finite-element-finite-time-difference method for solving the time-dependent coupled diffusion equations is firstly developed to serve as the forward model. Then, the inverse model for the simultaneous reconstruction of the fluorescent yield and lifetime images is proposed based on a self-normalized, time-binned data-type, in response to the temporal characteristics of the sate-of-the-art time-resolved measuring systems, and some issues associated to its implementation, such as the computation of the weighting matrix, linear inversion and regularization strategy, are described. The proposed algorithm is finally validated using simulated data for two-dimensional (2-D) scenarios, with its high performances in the improvement of image quality demonstrated in comparison with a featured-data scheme previously developed based on the generalized pulse spectrum technique (GPST)  and with the excitation-normalization.
2. Forward model - a finite-element-finite-difference solution to coupled time-domain diffusion equations
An extension of the numerical techniques for modeling the generic light propagation in turbid medium to the fluorescence regime is rather straight. Nevertheless, rephrasing the highlights of the method would be necessarily meaningful for the completeness and readability of the work since it thereafter lays the basis of the inverse model. Let the subscripts x and m denote the excitation and emission wavelengths, respectively. Then the propagation of both excitation and emission light in a heterogeneous, fluorescent turbid medium Ω, as externally excited by an ideal ultra-short point source (i.e. δ -shaped in space and time), can be modeled by the coupled time-domain diffusion equations as follows :
where ϕν(r,r s,t) (ν∈[x,m]) (t : time, r : position in Ω) is the temporally- and spatially-varying photon density corresponding to a δ -shaped source at position r s ; the optical properties involved are the absorption coefficient μaν(r), the reduced scattering coefficient μ′sν(r) and the diffusion coefficient κν(r,t)=c/[3μ′sν(r)] at the two wavelengths, respectively; the fluorescent parameters are the yield η(r) (product of the quantum efficiency and the absorption coefficient of the fluorophore) and the lifetime τ(r); e(τ,t)=e -t/τ U(t) with U(t) being a unit-step function; the operator ⊗ denotes the temporal convolution. We employ uniquely the Robin boundary condition for the above coupled equations. The measurable temporally-varying flux density, i.e., the temporal point spread function (TPSF), at the detection site ξd(d=1,2,…,D) on the boundary ∂Ω for the source site ζs(s=1,2,…,S), can be calculated by the Fick’s law under the Robin boundary condition 
where Rf≈0.53 is the internal reflection coefficient for air-tissue interface . For finite-element-method (FEM) solution to Eq.(1), ϕν(r,t) (r s is omitted for a brevity) is firstly approximated by a piecewise-polynomial function
where u(r)=[u 1(r),…,uN(r)]T and Φ ν(t)=[Φν(1,t),…,Φν(N,t)]T are vectors representing the shape functions and temporally-varying photon density at the N nodes of the FEM mesh. Then, a temporally-recursion matrix equation can be obtained below by employing the standard Galerkin-FEM procedure following by a full implicit time-domain-finite-difference (TDFD) scheme [17, 21] for a sequence of time points with an time interval of Δtf : Φ ν(i)=Φ ν(iΔtf), i=-1,0,1,…
where A ν, B and C are matrices of N×N with the same expression as given in , for both the excitation and emission wavelengths, i.e.
But Q ν(i) differs in form
where Φ′x(i)=[Φ′x(1,i),…,Φ′x(N,i)]T with , η(n) and τ(n) being the fluorescent yield and lifetime at the n-th node, respectively, and E(n,i)=e(τ(n),iΔtf); δ(i) is a Kronecker δ -function. Since the matrix (A ν+B)+C/Δtf is constant throughout the TDFD recursion in Eq. (4), the Choleski decomposition is needed only once during the solution process .
Without loss of the generality, we consider throughout this study a circular domain with a radius of R=15 mm and optical and fluorescent properties of μ(B) ax,m=0.035 mm-1 μ′(B) sx,m=1.0 mm-1, η(B)=0.001 mm-1 and τ(B)=500 ps ((B) denotes the background). The geometry and parameters are chosen for the potential applications of FMT in small animals as well as the status of the currently favorable red-shifted or NIR fluorescent probes [3, 8, 22, 23]. As an example of the proposed FEM-TDFD forward model, we assume that a small fluorescent region with a radius of r=3 mm and the same optical properties as the background is embedded at the center. The domain is excited by a δ- shaped laser pulse at a polar angle of θ=0° and the excitation re-emission and fluorescent emission are temporally measured at θ=180°, respectively. The FEM mesh contains 2400 triangle elements that join at 1261 nodes. Figure 1(a) shows the geometrical setup and the FEM mesh. In Case 1, the fluorescent inclusion has a fixed yield of η=0.002 mm-1 but different lifetimes of τ=560, 1000, 2800 ps, respectively, covering the range of the commonly-used fluorophore agents, such as Indocyanine Green (ICG) and Cy-series (Cy5.5, Cy3-B) etc. . Figure 1(b) shows the normalized TPSFs of the fluorescent emissions and the excitation re-emission at the detector site, calculated by the FEM-TDFD scheme. It is seen from these TSPFs last approximately from 2 ns for the excitation to 6 ns for the fluorescent lifetime of 2.8 ns. The calculation presents a rough hint that, for the efficient acquisition of the information for small-animal imaging, the measuring system might need to possess a temporal resolution less than 200 ps and a time range around 10 ns, though the detailed specifications is highly subjective to the practical situations. In Case 2, we have checked the influence of the TDFD-interval on the accuracy of the forward calculation. Figure 1(c) presents the TPSFs of the fluorescent emission calculated with a fixed fluorescent lifetime of 2.8 ns and different TDFD-intervals: Δtf=10, 20, 40, 80 ps. The results show that the relative error of the maximum TPSF points for the two successive steps of the TDFD-interval: Δtf=10 ps and 20 ps, is smaller than 10-3. It is requested that, for assuring an adequate accuracy of the forward calculation while avoiding an excess computational cost, a critical TDFD-interval, e.g., 20 ps for the case, should be suitably decided.
3. Inverse model - a self-normalized, full time-resolved scheme for simultaneous reconstruction of fluorescent yield and lifetime
An original implementation of the full time-resolved scheme uses all temporal samples that are necessarily generated during the time-resolved measurements. Let F m[p(r)]=[Fm(ξ1,ζ1,t 1,p(r)),…,Fm(ξD,ζS,tI,p(r))]T symbolize the forward operator that maps the image space p(r)=[η(r), τ(r)]T into the data space Γ m=[Γm(ξ1,ζ1,t 1),…,Γm(ξD,ζS,tI)]T, with D, S and I being the numbers of the sources, detectors and temporal samples, respectively, and subscript m emphasizing the emission wavelength. Due to the nonlinear dependency of the measurements on τ(r), the inverse model for the simultaneous recovery of both the yield and lifetime can be in general addressed by the Newton-Raphson scheme, based on the Taylor series expansion of F m[p(r)]. This forms the outer-loop of an iterative procedure
where J(p) is the Fréchet derivatives of F m[p(r)] with respect to η(r) and τ(r). Its discrete formula, i.e., the Jacobian matrix, can be efficiently calculated from a perturbation approach to the emission diffusion equation, firstly leading to an integral expression as follows, associating the parameter perturbations, δη(r) and δτ(r) with the incurred measurement change δΓ m,
where f η(τ,t)=e(τ,t)/τ and f τ(τ,η,t)=e(τ,t)(t-τ)η/τ3 ; gm(r,r′,t) is the Green’s function associated with the emission diffusion equation. Similarly, the fluorescent yield and lifetime distributions can be decomposed over the shape function basis of the finite elements: δη(r)=δη T u(r) and δτ(r)=δτ T u(r), where δη=[δη(1),δη(2),…,δη(N)]T and δτ=[δτ(1),δτ(2),…,δτ(N)]T with δη(n) and δτ(n) being the yield and lifetime perturbations at the n-th node, respectively. Then, the n-th entries of the η- and τ- regarding Jacobian matrices for source site ζd, detector site ξs and time ti, J η(ξd,ζs,ti) and J τ(ξd,ζs,ti), are given, respectively, upon the spatial discretization
where Ωn→j denotes the j-th element associated with the n-th node (i.e., the element having the n-th node as its common vertex), V(Ωn→j) is the volume (or area) of element Ωn→j, and , represents the mean value over element Ωn→j, which is approximated by the average of the values at the nodes of Ωn→j. The factor is dependent on the shape of the used finite element as well as the order of the shape function . As a linear shape function is employed W=3 and 4 for a triangle and a tetrahedron, respectively. The twofold temporal convolutions involved in the above expressions are carried-out in the Fast Fourier Transform that is more time-saving but storage-costly than the direct linear convolution.
Finally, the algebraic reconstruction technique (ART) is employed for a solution to Eq. (7), whose row-fashioned performance offers a nearly memory-independent and robust inversion for the large-scale and ill-posed linear system and forms the inner-loop of the iterative procedure . To prevent the over-reconstruction adversity with using the ART inversion, the inner-loop iterations are performed to numerate all the rows of the matrix equation (one row representing the measurement of a certain source-detector combination) only once. The regularization strategy in the ART-embedded algorithm is accomplished by truncating the outer-loop iterations with a criterion that the relative error norm between the measured and predicted samples stops decreasing for two successive updating stages . For the artifact reduction, a median filtering is applied to the adjacent nodes (i.e., those that belong to the same element as the node of interest) at the end of each out-loop iteration, as it has been done in our previous studies [9, 17, 25].
3.1 Time-bin νersion
In realistic regime, measurement of light signal essentially represents a discrete number of photons in a given time-interval, i.e., the so-called time-bin, which is determined either by the bandwidth of the instrumental impulse response in an analog mode, or by the jitter performance of the output pulses in a digital one, such as the transit-time-spread (TTS) of the PMT in the TCSPC detection, or the gate-width in a boxcar mode, such as a gated ICCD camera, and accounts for the physical temporal resolution of the time-resolved measuring systems . In this regard, it is practically meaningful to adopt statistically averaged samples within the time-bin as a data-type for incorporating the temporal response of the detectors. With a bin-width of T, the data-type is written as
where tk=T[rand(k)-0.5] with rand(k) being a random variable uniformly distributed in [0,1] and K is an associate number with the photon counting. Accordingly, the Jacobian matrices with regard to the time-bin data-type become
It should be pointed out that Eq. (10) only describes a realistic model of the data formation while the experimental data-set that are actually put into the image reconstruction process is no other than the temporally binned flux samples m(ζd,ξs,ti) spaced at an sampling-interval of Δtr, i.e., ti=iΔtr,i=I 1,I 1+1, ⋯,I 2, where integers I 1 and I 2 usually vary with the sourcedetector pairs to cover the most effective range of the individual TPSFs, and Δtr, is in general determined by either the electronic resolution of a TCSPC system or the delay-interval of a gated ICCD camera.
Use of intensity-related data-type such as the above temporally binned flux always incurs the requirement for scaling or eliminating the coupling coefficients of the detection system that link the realistic measurements to the model predictions. Normally, this can be done with three kinds of strategies: the auxiliary-parameter optimization that solves for unknown coupling coefficients as part of the image reconstruction problem [27–29], the so-called difference imaging scheme that cancels the coupling coefficients using the ratio of the emission measured after fluorescence inclusion (task) to the background one [9, 12] and the normalized Born approximation [4, 19] that is essentially based on the intensity-independent data-type, referred to as the Born ratio.
The strategy of auxiliary-parameter optimization has been explored mainly in simulated DOT scenarios and shown to be a promising way of correcting the coupling effects [27–29]. Nevertheless, since introducing coupling coefficients as auxiliary parameters during image reconstruction increases somewhat the uncertainty in the inverse process, the image quality might be adversely affected. Although the severity can be alleviated considerably with pre-conditioning, the effectiveness of this strategy in fluorescence imaging is still in need of a full investigation, especially for realistic situations.
The difference imaging scheme has also been intensively studied in DOT regime and also enables the elimination of the most schematic errors common to both the background measurement and the task one [25, 41–43]. However, this is often an unrealistic technique in the fluorescence imaging since measurement at emission wavelength in absence of fluorescent agent may not provide a signal with adequately high SNR . With regard to the above situations, the Born ratio method has been widely accepted as an effective means of eliminating the system calibration.
There are two ways of formulating the Born ratio. The first one is to normalize the measured emission flux to the excitation one as the data-type , with the assumption that either the scaling-factor (system gain) is constant over the emission and excitation wavelengths, or its spectroscopic difference well calibrated, i.e.,
where . In this case, Eq. (6) is simply reformulated as
where F̅m(ξd,ζs,ti,η k,τ k) is the model-calculation of the temporally binned flux at the emission wavelength; F (E) x(ξd,ζs) is the model-calculation of the integrated intensity at the excitation wavelength; Θx and Θm are the scaling factors of the system at the excitation and emission wavelengths, respectively. The more effective way is to construct a self-normalization that normalizes the temporally binned flux m(ξd,ζs,t) to its integrated intensity:
With this data-type the requirement for the scaling-factor calibration at the excitation and emission wavelengths is released. Accordingly, we have the counterpart of Eq. (13) as follows
where F (E) m(ξd,ζs,η k) is the model-calculations of the integrated intensity at the emission wavelength, respectively; is the η-regarding Jacobian Matrix for the emission intensity with its entry for the n-th node given by [12, 13]
where Em(ξd,r) is the emission Green’s function and Φx(r,ζs) the excitation photon density.
In addition to the scaling-factor cancellation, the other practical benefits of employing the Born ratio might include the allowance of a 2-D reconstruction scheme and the substantial enhancement of the algorithm robustness to the background heterogeneity, as have been well demonstrated in both DOT and FMT investigations [5, 30, 31]. The cost incurred, however, would be degradation in the probing sensitivity and noise-robustness .
In the following section, the above full time-resolved methodology will be fairly contrasted with a GPST-based algorithm that is developed within the similar framework , for the evaluation of its performances.
4. Methodology evaluations
We employ the same 2-D domain as in Section 2 for the evaluations. As shown in Fig. 2, it is assumed that 16 optodes with coaxially structured source- and detection- fibers are equally spaced around the boundary, whose output ends are coupled into either a 16 TCSPC channels with a electronic temporal resolution of ≤25 ps and a physical temporal resolution of <200 ps, or a gated ICCD device with a gate-width of ≥500 ps and a delay-interval adjustable from several picoseconds up to several nanoseconds. This assumption is in agreement with the specifications of the conventionally available TCSPC systems [16, 26, 32] and time-gated ICCD cameras [2, 14, 15], meaning that a bin-width of T=180 ps for the TCSPC or T≥500 ps for the gated ICCD would be reasonable. To facilitate the quantification, two circular targets with the same optical properties as those of the background are aligned along the Y-axis, symmetric to the origin. They have different radius, r1 and r2, as well as fluorescent parameter sets, (η1,τ1) and (η2,τ2), respectively. For any reconstructed images, their Y-profiles, i.e., the cut lines along the Y-axis with notations (y) and (y) for the yield and lifetime, respectively, are extracted to evaluate the measures. We firstly investigate the intrinsic performance of the methodology with noiseless data, and then check its noise-robustness using a noise model in accordance with the physics of the fluorescence generation, propagation and detection. To validate the robustness of the method on the numerical error of the forward model, a fine mesh with 3750 elements and 1951 nodes is used in the data-generator, as compared with a coarse mesh for the inversion-solver that contains 2400 elements and 1261 nodes. Also, to assure an adequately high signal-to-noise ratio (SNR), we use the temporal samples ranging from 10% of the maximum intensity on the rising edge to 10% of the maximum intensity on the falling edge for each of the source-detector pairs. It is worth notice that some researchers have addressed the optimization of temporal measurements using a singular-value-decomposition analysis and shown a redundancy-reduction mode by obtaining an optimal result from the rise and the peak portions of the TPSFs .
We firstly introduce the spatial-resolution measure R χ(χ∈[ , ]) as the evaluation measure
Then the yield (η)- and lifetime (τ)-images are simultaneously reconstructed and the measures calculated from their Y-profiles for a varying center-to-center separation (CCS), using both the full time-resolved and GPST schemes. The two circular targets have the same radius of r 1=r 2=2 mm and fluorescent contrast to the background of 5:1 and 3:1 for the yield and lifetime, respectively, i.e., η1=η2=0.005 mm-1 and τ1=τ2=1500 ps. To investigate the effects of temporal features of the detectors on the reconstruction, a few of different time-bin widths and sampling intervals are used in the reconstruction. Figure 3 shows the e R χ as functions of the CCS for different sampling-intervals and different bin-widths, respectively. To better demonstrate the superiority of the time-resolved method to the featured-data one, the results obtained with the GPST algorithm are presented for a comparison. Figure 4 illustrates the reconstructed images using the proposed scheme for a fixed sampling-interval of tr=20 ps and a series of the bin-widths: T=20, 180, 500 and 1000 ps, while Figure 5 shows the reconstructed η - and τ -images using the proposed method for a fixed bin-width of T=180 ps and different sampling-intervals of tr=20, 100, and 180 ps, both at a fixed CCS=11 mm. It is noted from the results that the sampling-interval greatly affects spatial resolution of the reconstruction, and at a small sampling, e.g., tf=20 ps, the moderate increase of the bin-width (<500 ps) has no evidently adverse influence on the image quality.
We simply calculate the reconstructed fluorescent contrast Q χ for the evaluation
Similarly, we quantitatively investigate the performance by simultaneously reconstructing the η - and τ -images of the same scenario as the previous case using both the schemes for a fixed target CCS of 15 mm and a range of target fluorescent contrasts to the background from 2:1 to 8:1 for both the yield and lifetime, i.e., (ηi,τi)=(j0.001 mm-1, j500 ps), i=1,2; j=2, 3,…,8. Figure 6 calculates the reconstructed fluorescent contrast Q χ as a function of the target fluorescent contrast (i.e., the ratio of the maximum to the minimum in the target Y-profile) for different sampling-intervals and different bin-widths, respectively. For contrasting between the two schemes, the results from the GPST method are also included. The similar conclusion can be obtained from the results: the yield and (low-contrasted) lifetime reconstructions with high quantitativeness can be attained at a small sampling interval of tf=20 ps, and are insensitive to an increase in the bin-width up to 500 ps.
The size-contrast performance S χ as a measure of the fidelity of the methodology to distinguish the difference in the target size between the two targets, is defined as
where and are the reconstructed maximum values of the target 1 and 2, respectively; χ(min)=min[χ(y)] is the minimum value of the Y-profile; y (l) i and y (u) i (i=1,2) are the y-coordinates of the half-maximum of the target 1 and 2, respectively. Again, the same scenario as before is used with a fixed d CCS=15 mm for the two circular targets and fluorescent contrast of 5:1 and 3:1 for the target yield and lifetime, respectively. Figure 7 gives the S χ as a function of the target size-contrast (defined as the ratio of the diameter of the target 1 to that of the target 2), ranging from 1 to 3, i.e., r 1,2=2±0.25 j mm, j=0,1,…,5. Similarly, the calculations are made for different sampling-intervals and different bin-widths and compared to the GPST results. Again, we have the similar observation as the above that the measure is adversely affected by the increase in the sampling-interval and nearly unchanged with the bin-width under 500 ps.
The measure, G χ, is used to quantify the ability for the method to reconstruct the gray difference
where χ(max) 1, χ(max) 2 and χ(min) are afore-defined. We quantify the performance by calculating G χ as a function of the target grayscale difference at a fixed baseline of η=0.005 mm-1 and τ=1500 ps, i.e., the mean fluorescent parameters of the two targets with the radius of r 1=r 2=2 mm and CCS=15 mm, as shown in Fig. 8. In addition, we also evaluate G χ as a function of the target η-baseline at a fixed target lifetime baseline of τ=1500 ps, a target grayscale difference of 30% and a target radius of r 1=r 2=2 mm, and as a function of the target radius at a fixed target baseline of η=0.005 mm-1 and τ=1500 ps and a target grayscale difference of 30% (the results are not given for the brevity). The evaluations are contrasted for different sampling-intervals and different bin-widths, respectively, as well as compared with the GPST results. For better demonstrating the effect of the sampling-interval selection on the reconstruction fidelity of the grayscale difference, Figure 9 illustrates the reconstructed images of the above first scenario with tr=20, 100, and 180 ps, respectively, at a fixed bin width of T=180 ps. It is seen that from Fig. 8, for the yield reconstruction an optimal result is achieved at a sampling-interval of tf=20 ps, as compared to the GPST method that over-estimates the measure, and to the proposed one with a larger time-interval that under-estimates the true value. As to the lifetime reconstruction, the GPST method gives an over-estimation while the proposed one underestimation with the worst result generated at tf=20 ps. The observation is further manifested in Fig. 9.
Generally, two types of noise sources are dominant in the measured signals. The first is due to thermal-emission in the detector electronics. This type of noise is independent of the measured signal and well modeled as a sequence of independent identically distributed Gaussian random variables. The second one, known as shot-noise, originates from the quantum fluctuations of the detected photons, and is typically modeled as a Poisson distributed variation in the measured flux density. This model can be further approximated with independent, but not identically distributed Gaussian random variables, with the SNR equal to the square root of the (averaging) number of the detected photons at time t in a given measurement (integration) time . With an expectation that measured flux density will be reasonably high to exceed the thermal noise floor of the detectors, we will be able to employ the following noise model for characterizing the noise-robustness of the full time-resolved scheme
where Npm is the number of the detected photons corresponding to min[Γm(ξd,ζs,t)]. We introduce the mean value of SNR(ξd,ζs,t) over the time range, to depict the overall noise-level in the measured TPSFs.
Figure 10 presents the simultaneously reconstructed η - and τ -images at a mean SNR level of 25 dB, for different bin-widths of T=20, 180, and 500 ps, respectively, at tr=20 ps, where a standard target settings with η - and τ -contrasts of 5:1 and 3:1, respectively, a radius of r 1=r 2=2 mm and a CCS of 15 mm are used. An overall improvement in the image quality can be evidently observed at T=500 ps
Contrasting to the featured-data method, e.g., the GPST scheme, perhaps, the biggest gain with using full time-resolved data lies on the significant improvement in the spatial resolution of both the yield and lifetime images. This is achieved for both the yield and lifetime reconstructions at a sufficiently narrow sampling-interval, e.g., tr=20 ps, as clearly seen in Fig. 3. The observations are further enforced with the reconstructed images in Fig. 4 and Fig. 5, where the resultant images exhibit a quality consistency for T≤500 ps and an obvious quality degradation with the increasing sampling-interval. According to Fig. 6, the proposed full time-resolved scheme significantly improves quantitativeness of the yield reconstruction but evidently degrades the quantitation in the lifetime image for a target-to-background contrast higher than 4:1. The latter disadvantage is primarily ascribed to the increasingly saturating sensitivity of the detected signals to the long-lifetime fluorescent variation. As illustrated in Fig. 7, the size-contrast reconstructions with the proposed method is slightly under-estimated for the yield reconstruction and slightly over-estimated for the lifetime one, whereas the GPST method presents an severely over-estimated reconstruction for the yield and a slight over-estimation for the lifetime, in terms of the measure S χ. Similarly, an optimal reconstruction fidelity of the target size contrast is found at a small time-interval of tr=20 ps. The evaluations on the grayscale performance of the method show that the reconstruction with the full time-resolved method is under-estimated for both the yield and lifetime, contrary to those with the GPST, as shown in Fig. 8. Further evaluations on this performance as functions of the target contrast and size, respectively, show that with the increase in the yield contrast and target radius, the grayscale reconstruction of the yield gets successively improved while the reconstructed grayscale of the lifetime degrades slightly. Roughly, the results conform to our previous observations in DOT . Beside the improvement in image quality with the full time-resolved scheme, another important observation from the above evaluations is that the quality of the image reconstruction is largely dependent on the sampling-interval but the bin-width, and the superior performances of the full time-resolved scheme over the featured-data one is achievable as the sampling-interval gets sufficiently small, again as demonstrated in Fig. 9. From an information point of view, it is very natural for the performance to be interpreted because a wide bin with a small sampling-interval does not mean any loss of the information content embedded in the time-resolved data. In practice, however, this feature is significantly beneficial to the choice of a time-resolved detection mode, in particular, allowing for the use of a state of the art ICCD camera with a gating-time of several hundred picoseconds for a densely spatial sampling, or the channel-overlay mode in a TCSPC system for a SNR improvement of the data-set. As demonstrated in Fig. 10, properly broadening the bin-width of detecting system at a sufficiently small sampling-interval leads to an effective enhancement in the noise-robustness of the reconstruction, while still reserves the advantages of the full time-resolved scheme.
The method described considers the samples with the uniform optical absorption and scattering properties. This situation is never observed in practice. Therefore, developing robust methods that work at highly heterogeneous background is realistically important. Soubret et. al. have investigated the influence of local variations in the optical properties on the accuracy of the fluorescence reconstruction and demonstrated that, the reasonable quantification is relatively insensitive to a degree of background heterogeneities, with the excitation-normalized Born approach . As for the self-normalized scheme proposed here, despite some degradation in the performances compared with the excitation-normalization, the relative independency of the normalized methods on the local variations of the background optical properties is still retained upon the reasonable assumption that the optical heterogeneities are prone to evidently affecting the absolute intensity of the time-resolved fluorescent signals but its relative variations . Nevertheless, a hybrid modality of FMT and DOT that is capable of providing an inversion scheme in terms of realistic optical anatomy is greatly desirable for the enhancement of the reconstruction sensitivity and accuracy, as well as for the universal applicability of the proposed methodology.
It is noted here that, for the calculation simplicity and methodology demonstration, the data generation and image reconstruction in this study are all based on the photon diffusion theory that has been well established in the field of the diffuse light imaging [21, 25]. For the validation of the method robustness on numerical error of the forward calculations, a finer FEM mesh is used for the generation of the simulated data, for which a noisy data model is further developed. In coming investigations, it is requested that more convincing validations be carried on with more realistic data-sets, such as those from either Monte-Carlo simulations or experiments, to alleviate the “inverse crime” issue that may arise from using the identical mathematical model in the forward and inverse calculations. Inevitably, a further degradation in image quality would be observed due to the model mismatch between the mathematics of the inversion strategy and the physics of the measurement formation.
This study mainly considers the Poisson-distributed shot noise, which is further approximated with a photon-level dependent Gaussian stochastic process, as afore-mentioned . Although it is somewhat empirical, the model is practically reasonable for the general evaluation of the method noise-robustness. For working with the full time-resolved data, a systematic error of measurement that might markedly influence the reconstruction originates from the temporal offsets of the data-sets, as shown in Fig. 11, where a simultaneous reconstruction is performed with noiseless data for the aforementioned background embedding at (x=0 mm, y=5.5 mm) a circular target with a radius of r=3 mm and η - and τ - contrasts of 5:1 and 3:1, respectively. According to the results a margin of the time origin within ±80 ps is required for a reasonable output, which is somewhat rigorous in practice. Normally, the solution to this issue is strongly dependent on the availability of the calibrating tools that directly measure the absolute time-origins of the system channels with an acceptable accuracy. Although some effective efforts have been made for this direct strategy , the solution remains less satisfactory due to either the imperfectness of the methodology or the instability of the system, in particular, for the absolute imaging cases, such as FMT. Alternatively, an indirect method, which avoids the accurate determination of the absolute time-origin by firstly convoluting the model predictions with the measured instrumental response functions and then matching the calibrated data-sets to the measured ones, would be promising. A further evaluation indicates that increasing the bin-width will enhance the robustness of the reconstruction to the temporal offsets to some extend.
Figure 12 illustrates the Y-profiles of the reconstructed η- and τ - images on a scenario similar to in Fig. 4 but with a target CCS of 10 mm, for a brief comparison between the excitation- and emission-normalized formulations. The results with the emission-normalization show slightly inferior η-resolution and superior τ -resolution to those obtained with the excitation-normalization, whereas the reconstructed quantitativeness is almost equivalent for both the schemes. This observation is significantly meaningful in practice with considering the benefits that the self-normalized version potentially offers.
Normally, the spatial resolution and reconstruction quantitativeness are two crucial performances of imaging modalities. However, due to the high illposedness presented in the model-based diffuse optical imaging techniques, the inversion process tends to induce crosstalks between the target size and contrast, i.e., the difference in size is often converted in part to difference in contrast, and vice versa. Therefore, the size-contrast and gray-resolution measures introduced in this paper are necessarily required for the completeness and objectiveness of the methodology evaluations.
We have proposed a general framework for the TD fluorescence-DOT, aiming at simultaneous reconstruction of the fluorescent yield and lifetime images using the full time-resolved data that is self-normalized and time-binned for simplifying implementation and adapting to the system reality. The methodology was validated with simulated data for its performances in the spatial-resolution, quantitative accuracy, reconstruction fidelity of size and grayscale, as well as the noise-robustness. The results have showed an overall superiority of the full time-resolved scheme to the featured-data one, as expected. In addition, two important observations are concluded:
- Compared to the Born formulation where the ratio of emission data to excitation one is employed for the elimination of the system scaling, the self-normalization scheme possesses additional advantages of enhancing the robustness to system fluctuations and exempting spectroscopic characterization of the system, whereas the performance degradation incurred due the loss in the intensity information is unconspicuous;
- A noise-robust image reconstruction might be attained with reasonable fidelity by properly increasing the bin-width, e.g., up to 500 ps, while remaining the sampling-interval enough narrow, e.g., down to 25 ps. This feature technically allows for the application of a state of the art ICCD for dense spatial sampling, and offers a favorable way of improving the SNR.
In summary, a time-domain mode is currently trending towards a focus of the future technological development, with rapid advances in ultrafast detection hardware and pulsed laser diodes. The time-resolved data obtained with the mode, however, needs to be properly organized to fully dig-up the embedded information for a high-fidelity reconstruction as well as for a balance among computation cost, noise-robustness and practical feasibility. The methodology described in this paper offers a flexible approach to the goal. The experimental validation of the proposed scheme is on-going now and the preliminary results will be reported on later.
Feng Gao and Huijuan Zhao acknowledge the funding supports from the National Natural Science Foundation of China (60578008, 60678049), the National Basic Research Program of China (2006CB705700) and Tianjin municipal government of China (07JCYBJC06600). Yukio Yamada thanks the support from Grant-in-Aid for Scientific Research (B) (19360098) by Japan Society for the Promotion of Science (JSPS).
References and Links
4. V. Ntziachiristos, J. Ripoll, L. H. V. Wang, and R. Weissleder, “Looking and listening to light: the evolution of whole-body photonic imaging,” Nat. Biotech. 23, 313–320 (2005). [CrossRef]
6. X. Cong and G. Wang, “A finite-element-based reconstruction method for 3D fluorescence tomography,” Opt. Express 139847–9857 (2005), http://www.opticsinfobase.org/abstract.cfm?URI=oe-13-24-9847. [CrossRef] [PubMed]
7. A. T. N. Kumar, S. B. Raymond, G. Boverman, D. A. Boas, and B. J. Bacskai, “Time-resolved fluorescence tomography of turbid media based on lifetime contrast,” Opt. Express 14, 12255–12270 (2006), http://www.opticsinfobase.org/abstract.cfm?URI=oe-14-25-12255. [CrossRef] [PubMed]
8. S. Lam, F. Lesage, and X. Intes X, “Time-domain fluorescent diffuse optical tomography: analytical expressions,” Opt. Express 13, 2263–2275 (2005), http://www.opticsinfobase.org/abstract.cfm?URI=oe-13-7-2263. [CrossRef] [PubMed]
9. F. Gao, H. J. Zhao, Y. Tanikawa, and Y. Yamada, “A linear, featured-data scheme for image reconstruction in time-domain fluorescence molecular tomography,” Opt. Express 14, 7109–7124 (2006), http://www.opticsinfobase.org/abstract.cfm?URI=oe-14-16-7109. [CrossRef] [PubMed]
10. H. Jiang, “Frequency-domain fluorescent diffusion tomography: a finite-element-based algorithm and simulation,” Appl. Opt. 37, 5337–5343 (1998). [CrossRef]
11. A. D. Klose, V. Ntziachristos, and A. D. Hielscher, “The inverse source problem based on the radiative transfer equation in optical molecular imaging,” J. Comput. Phys. 202, 323–345 (2005). [CrossRef]
12. J. Lee and E. M. Sevick-Muraca, “Three-dimensional fluorescence enhanced optical tomography using referenced frequency-domain photon migration measurements at emission and excitation wavelengths,” J. Opt. Soc. Am. A 19, 759–771 (2002). [CrossRef]
13. R. Roy, A. B. Thompson, A. Godavarty, and E. M. Sevick-Muraca, “Tomographic fluorescence imaging in tissue phantom: A novel reconstruction algorithm and imaging geometry,” IEEE Trans. Med. Imaging 24, 137–154 (2005). [CrossRef] [PubMed]
14. V. Y. Soloviev, K. B. Tahir, J. McGinty, D. S. Elson, M. A. A. Neil, P. M. W. French, and S. R. Arridge, “Fluorescence lifetime imaging by using time-gated data,” Appl. Opt. 46, 7384–7391 (2007). [CrossRef] [PubMed]
16. M. Brambilla, L. Spinelli, A. Pifferi, A. Torricelli, and R. Cubeddu, “Time-resolved scanning system for double reflectance and transmittance fluorescence imaging of diffusive media,” Rev. Sci. Instrum. 79, 013103 (2008). [CrossRef] [PubMed]
18. T. J. Farrell and M. S. Patterson, “Diffusion modeling of fluorescence in tissue,” in Handbook of Biomedical Fluorescence, Mycek MA and Pogue BW eds., Marcel Dekker, New York (2003). [CrossRef]
19. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15, R41–93 (1999). [CrossRef]
20. W. G. Egan and T. W. Hilgeman, Optical Properties of Inhomogeneous Materials (Academic Press, New York1979).
22. S. Achilefu, P. R. Dorshow, J. E. Bugaj, and R. Rajapopalan, “Novel receptor-targeted fluorescent contrast agents for in vivo tumor imaging,” Invest. Radiol. 35, 479–485 (2000). [CrossRef] [PubMed]
23. K. Licha, “Contrast agents for optical imaging,” Topics in Current Chemistry 222, 1–29 (2002). [CrossRef]
24. O. C. Zienkiewicz and R. L. Taylor, The Finite Element Methods Volume 1 5th ed., (Elsevier Pte Ltd., Singapore2004).
25. H. J. Zhao, F. Gao, Y. Tanikawa, and Y. Yamada, “Time-resolved diffuse optical tomography and its application to in vitro and in vivo imaging,” J. Biomed. Opt. 12, Art. No. 062107 (2007). [CrossRef] [PubMed]
26. W. Becker, Advanced time-correlated single photon counting techniques (Springer-Verlag, Berlin2005). [CrossRef]
27. D. A. Boas, T. Gaudette, and S. R. Arridge, “Simultaneous imaging and optode calibration with diffuse optical tomography,” Opt. Express 8, 253–270 (2001). [CrossRef]
28. S. Oh, A. B. Milstein, R. P. Millane, C. A. Bouman, and K. J. Webb, “Source-detector calibration in threedimensional Baysian optical diffusion tomography,” J. Opt. Soc. Am. A 19, 1983–1993 (2002). [CrossRef]
29. T. Tarvainen, V. Kolehmainen, M. Vauhkonen, A. Vanne, A. P. Gibson, M. Schweiger, S. R. Arridge, and A. P. Kaipio, “Computational calibration method for optical tomography,” Appl. Opt. 44, 1879–1888 (2005). [CrossRef] [PubMed]
30. J. C. Hebden, A. Gibson, T. Austin, R. Yusof, N. Everdell, D. T. Delpy, S. R. Arridge, J. H. Meek, and J. S. Wyatt, “Imaging changes in blood volume and oxygenation in the newborn infant brain using three-dimensional optical tomography,” Phys. Med. Biol. 49, 1117–1130 (2004). [CrossRef] [PubMed]
31. H. J. Zhao, F. Gao, Y. Tanikawa, K. Homma, and Y. Yamada, “Time-resolved optical tomographic imaging for the provision of both anatomical and functional information about biological tissue,” Appl. Opt. 43, 1905–1916 (2005). [CrossRef]
32. F. E. W. Schmidt, M. E. Fry, E. M. C. Hillman, J. C. Hebden, and D. T. Delpy, “A 32-channel time-resolved instrument for medical optical tomography,” Rev. Sci. Instrum. 71, 256–265 (2000). [CrossRef]
33. F. Gao, H. J. Zhao, Y. Tanikawa, K. Homma, and Y. Yamada, “Influences of target size and contrast on near infrared diffuse optical tomography - a comparison between featured-data and full time-resolved schemes,” Opt. Quantum Electron. 37, 1287–1304 (2005). [CrossRef]
34. A. Soubret, J. Ripoll, and V. Ntziachristos, “Accuracy of fluorescent tomography in the presence of heterogeneities: study of the normalized Born ratio,” IEEE Trans. Med. Imaging 24, 1377–1386 (2005). [CrossRef] [PubMed]
35. F. Gao, P. Poulet, and Y. Yamada, “Simultaneous mapping of absorption and scattering coefficients from full three-dimensional model of time-resolved optical tomography,” Appl. Opt. 39, 5898–5910 (2000). [CrossRef]
36. E. M. C. Hillman, J. C. Hebden, F. E. W. Schmidt, S. R. Arridge, M. Schweiger, H. Dehghani, and D. T. Delpy, “Calibration techniques and datatype extraction for time-resolved optical tomography,” Rev. Sci. Instrum. 71, 3415–3427 (2000). [CrossRef]