OSA's Digital Library

Optics Express

Optics Express

  • Editor: Michael Duncan
  • Vol. 14, Iss. 16 — Aug. 7, 2006
  • pp: 7109–7124
« Show journal navigation

A linear, featured-data scheme for image reconstruction in time-domain fluorescence molecular tomography

Feng Gao, Huijuan Zhao, Yukari Tanikawa, and Yukio Yamada  »View Author Affiliations


Optics Express, Vol. 14, Issue 16, pp. 7109-7124 (2006)
http://dx.doi.org/10.1364/OE.14.007109


View Full Text Article

Acrobat PDF (507 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

Fluorescence diffuse optical tomography (DOT) has attracted many attentions from the community of biomedical imaging, since it provides effective enhancement in imaging contrast. This modality is now rapidly evolving as a potential means of monitoring molecular events in small living organisms with help of molecule-specific contrast agents, referred to as fluorescence molecular tomography (FMT). FMT could greatly promote pathogenesis research, drug development, and therapeutic intervention. Although FMT in steady-state and frequency-domain modes have been heavily investigated, the extension to time-domain scheme is imminent for its several unique advantages over the others. By extending the previously developed generalized pulse spectrum technique for time-domain DOT, we propose a linear, featured-data image reconstruction algorithm for time-domain FMT that can simultaneously reconstruct both fluorescent yield and lifetime images of multiple fluorephores, and validate the methodology with simulated data.

© 2006 Optical Society of America

1. Introduction

2. Methodology

Analogous to DOT, the task of image reconstruction in FMT is also expressed as an inverse issue for a given photon-migration model. We describe in this section the FEM frameworks for both the forward and inverse models of time-domain FMT, aiming at simultaneous reconstruction of fluorescent yields and lifetimes of multiple components from the time-resolved measurements.

2.1 Forward Model

In the field of diffuse light imaging, the diffusion equation, i.e. the P 1 approximation to the radiative transfer equation, has been prevalently employed as the photon-migration model, which is mathematically tractable either analytically or numerically [18

18. S.R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15, R41–93 (1999). [CrossRef]

]. For time-domain FMT, where an ultra-short pulsed laser at the excitation wavelength is launched into tissue and the resulting time-dependent fluorescent emission on the boundary is measured by a time-resolved system, the forward model that governs the relationship between the excitation and emission light propagations in fluorescent turbid medium is thus given by the time-dependent coupled diffusion equations. To apply the GPST, the time-dependent forward model is converted into the complex frequency domain by using Laplace transform

{[Dx(r)μax(r)cp]Φx(r,rs,p)=δ(rrs)[Dm(r)μam(r)cp]Φm(r,rs,p)=Φx(r,rs,p)ημaf(r)[1+pτ(r)]
(1)

where subscripts x and m denote the excitation and emission wavelengths, respectively; Φv(r,r s,p) = 0+ Φv(r,r s,t)e -pt dt (v∈[x,m]) is the Laplace transform of the time-dependent photon density Φv(r,r s,t) with a complex transform-factor p ; the optical parameters involved are the absorption coefficient μav(r), the reduced scattering coefficient μ′sv(r) and the diffusion coefficient Dv (r,t) = c/[3μ′sv(r)]; the fluorescence parameters are the fluorescent yield ημ(r) and lifetime τ(r). These quantities, in general, are functions of the position vector r. We employ uniquely the Robin boundary condition for the above equations

Φv(r,rs,p)+2KD(r)n∙Φv(r,rs,p)rΩ=0
(2)

where K = (1 + Rf )/(1 - Rf ) and Rf ≈ -1.4399n -2 + 0.7099n -1 + 0.6681 + 0.0636n is the internal reflection coefficient at the air-tissue boundary, with n being the relative refractive index of tissue to air [28

28. W.G. Egan and T.W. Hilgeman, Optical Properties of Inhomogeneous Materials, (Academic, New York, 1979).

]. The measurable flux, i.e., the data-type at the boundary site ξd (d = 1,2,…,D)) and for the source site ζs (s = 1,2,…,S), can be calculated by the Fick’s law with consideration of Eq.(2)

Γv(ξd,ζs,p)=cΦv(ξd,ζs,p)(2K)
(3)

Eq. (1) can be efficiently solved with the Galerkin FEM by expanding Φv(r,r s,p) into finite elements: Φv(r,r s,p) ≈ Σn=1N Φv(n,p) un (r) = Φ v(p)T u(r) with u(r) = [u 1(r), u 2(r), ⋯, uN (r)T and Φ v(p) = [Φv(1,p),Φv(2,p),⋯,Φv(N,p)]T being the shape functions and the Laplace-transformed photon density at the N nodes of the FEM mesh, respectively, resulting in following matrix equation

(Av+B)Φv(p)=Qv
(4)

The elements of the coefficient matrices A v, B and Q v are given below

{Avij=Ω{Dv(r)ui(r)uj(r)+[μav(r)c+p]ui(r)uj(r)}drBij=cΩui(r)uj(r)dr(2K)Cij=Ωui(r)uj(r)drQvit={Ωui(r)δ(rrs)dΩ=ui(rs)v=xj=1NCijΦxjpημaf(j)[1+pτ(j)]v=m
(5)

where i and j index the N nodes of the mesh; ημaf(j) and τ(j) the fluorescent yield and lifetime at the j-th meshing node, respectively.

2.2 Inverse model

According to Eq. (1) and by applying the Fick’s Law at emission wavelength, we have an integral equation as follows

{Γm(ξd,ζs,p)=ΩGm(ξd,r,p)ΦxrζspxrpdΩxrp=ημaf(r)[1+pτ(r)]
(6)

where Γmd, ζs,p) is the Laplace transform of the transient emission flux measured at boundary site ξd and for excitation site ζs, and Gmd, r, p) the flux at ξd for a source at r, predicted by the diffusion equation at the emission wavelength

{[Dm(r)μam(r)cp]Φm(r,r',p)=δ(rr')Φm(r,r',p)+2KD(r)n∙Φm(r,r',p)rΩ=0Gm(ξd,r',p)=(cK2)Φm(r,r',p)r=ξd
(7)

Let x(r,p)≈ Σn=1N xn (p)un (r) = x T(p)u(r) , where x(p) = [x 1(p), x 2(p),…, xN (p)]T , Eq.(6) can be discretized into a matrix equation

Γ(p)=W(p)x(p)
(8)

where

Γ(p)=[Γm(ξ1,ζ1,p),Γm(ξ2,ζ1,p),,Γm(ξD,ζS,p)]T
(9)

and

W(p)=[W(ξ1,ζ1,p,1),W(ξ1,ζ1,p,2),,W(ξ1,ζ1,p,N)W(ξ2,ζ1,p,1),W(ξ2,ζ1,p,2),,W(ξ2,ζ1,p,N)W(ξD,ζS,p,1),W(ξD,ζS,p,2),,W(ξD,ζS,p,N)]
(10)

with the element given by

Wds(ξd,ζs,p,n)=ΩnḠm(Ωn)(ξd,p)Φ̄m(Ωn)(ζs,p)Ωnun(r)dΩ
(11)

For solving Eq. (8) that is in general of large-scale, underdetermined and ill-posed, a Kacmarcz method, commonly known as the Algebraic Reconstruction Technique (ART) might be very efficient [18

18. S.R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15, R41–93 (1999). [CrossRef]

, 22

22. F. Gao, P. Poulet, and Y. Yamada, “Simultaneous mapping of absorption and scattering coefficients from a three-dimensional model of time-resolved optical tomography,” App. Opt 39, 5898–5910 (2001). [CrossRef]

, 29

29. A.C. Kak and M. Slaney, Principle of Computerized Tomographic Imaging, (IEEE Press, New York, 1988).

]

{xk+1(p)=xk(p)+λΓ((kmodSD)+1)(p)W(kmodSD)+1)(p)xk(p)[W(kmodSD+1)(p)][W(kmodSD+1)(p)Tk=0,1,2,(MSD1)
(12)

where Γ (i)(p) is the i-th element of Γ(p) and W (i)(p) the i-th row of W(p), with (i mod j) denoting the modulus after i divided by j ; M is referred to as the ART-circulation number. The initial point x 0(p) is set to what reflects the morphology of the background fluorescence as closely as possible. In principle, the ART sequentially projects a solution estimate onto the hyperplanes defined by the individual rows of the linear system. The relaxation parameter λ has a significant influence on the reconstruction and has been proven in a range of [0, 2] to make the algorithm converge to a point on the intersection of the governing equations that is nearest to the initial point [29

29. A.C. Kak and M. Slaney, Principle of Computerized Tomographic Imaging, (IEEE Press, New York, 1988).

]. The regularization strategy in the ART-based algorithms is accomplished by limiting the number of the iterations, whose choice is task-dependent and mandatory in presence of noise. The primary advantage of this implicit inversion method over the other schemes for solving underdetermined linear systems, such as the Levenberg-Marquardt, the truncated Singular Value Decomposition and the Conjugate-Gradient algorithms, is its near independence of memory occupation, since the rows of the weight matrix are successively employed during the solution process. However, the accuracy of the ART strongly relies on the initial guess of the unknowns. Sometimes, the prior information on the background as well as the targets is necessarily required to attain an image reconstruction of high-quality. To suppress the artifacts arising in the above ART inversion, a median filter operating on the adjacent nodes (i.e., those that belong to the same element as the output node) is employed at the end of each circulation of the ART inversion.

2.2.1 One-component case

In the normal one-component case, one kind of biochemical molecule of interest is targeted by its specific probe, and a set of two unknown distributions (one component), i.e. the fluorescent yield ημaf(r) and lifetime τ(r) of the probe, is to be recovered. This can be explicitly done from the images of x(r,p 1) and x(r,p 2) by employing a pair of transform-factors: p 1 and p 2, in the Laplace transforms

{ημaf(r)=(p1p2)x(r,p1)x(r,p2)[p1x(r,p1)p2x(r,p2)]τ(r)=[x(r,p1)x(r,p2)][p1x(r,p1)p2x(r,p2)]
(13)

2.2.2 Multi-component case

xrp=i=1Ncημafi(r)[1+pτi(r)]
(14)

where Nc is the number of the components under investigation. The recovery of the 2Nc unknown distributions in Eq. (14) can in principle be achieved by using at least 2Nc transform-factors in the Laplace transforms: p 1, p 2,…, p2Nc . For the two-component case, this leads to a set of four joint equations as follows

{ημa1(r)[1+p1τ1(r)]+ημa2(r)[1+p1τ2(r)]=x(r,p1)ημa1(r)[1+p2τ1(r)]+ημa2(r)[1+p2τ2(r)]=x(r,p2)ημa1(r)[1+p2Ncτ1(r)]+ημa2(r)[1+p2Ncτ2(r)]=x(r,p2Nc)
(15)

where x(r,pi ) (i = 1,2,…,2Nc ) is obtained from solving Eq. (8) for each of the transform-factors. Let

{α1(r)=ημaf1(r)+ημaf2(r)α2(r)=ημaf1(r)τ2(r)+ημaf2(r)τ1(r)α3(r)=τ1(r)+τ2(r)α4(r)=τ1(r)+τ2(r)
(16)

, Eq. (15) can be transformed into the following matrix equation

x(r)=M(r)α(r)
(17)

with

M(r)=[1,p1,p1x(r,p1),p12x(r,p1)1,p2,p2x(r,p2),p22x(r,p2)1,p2Nc,p2Ncx(r,p2Nc)p2Nc2x(r,p2Nc)]
(18)
α(r)=[α1(r),α2(r),α3(r),α4(r)]T
(19)
x(r)=[x(r,p1),x(r,p2),,x(r,p2Nc)]T
(20)

By solving Eq. (16) for α(r), we can further obtain the four unknown distributions of the fluorescence parameters, i.e., ημaf1(r), τ1(r), ημaf2(r) and τ2(r), in terms of Eq. (16)

{τ1(r)=[α3(r)α32(r)4α4(r)]2τ2(r)=[α3(r)+α32(r)4α4(r)]2ημaf1(r)=[α1(r)τ1(r)α2(r)][τ1(r)τ2(r)]ημaf2(r)=[α1(r)τ2(r)α2(r)][τ1(r)τ2(r)]
(21)

For the case of Nc > 2, a similar procedure to the two-component case can be followed, which is mathematically more complex. It should be noted that within the above FEM framework Eq. (17) and Eq. (21) need to be solved successively for each node of the mesh, given the values of x(r) at the positions of the discrete nodes.

Although, according to Eq. (17), 2Nc , real transform-factors are theoretically enough for the reconstruction of the 2Nc , independent unknowns in a Nc -component case, this may bring about considerable errors in the images due to the evident underestimation of x(r, p) in Eq. (8). To overcome the adversity, more than 2Nc real transform-factors are in general employed, making Eq. (17) both overestimated and ill-posed. We therefore introduce an extended Levenberg-Marquardt method for the reliable solution, which is expressed as an optimization problem with the Tikhonov-Miller regularization as the follows [30

30. F. Gao, H. Niu, H. Zhao, and H. Zhang, “The forward and inverse models in time-resolved optical tomography imaging and their finite-element method solutions,” Image and Vision Computing 16, 703–712 (1998). [CrossRef]

]

α(r)=argmin{x(r)M(r)α(r)2+βR(r)[α(r)αB(r)]2}
(22)

where α B(r) is the value of α(r) calculated for the known background fluorescent properties at each node. The regularization matrix R(r)is chosen in such a way that R T R (RR T) and M TM (MM T) commute each other, meaning that the two matrices have the same complete set of eigenvectors and the reversely-ordered eigenvalues, i.e., M=μkukvkTandR=k=1RMμRMkμkvkT, where RM is the rank of M, u k and v k are the eigenvectors of MM T and M T M belonging to the nonzero eigenvalue μk. The solution to the above optimization problem can be readily found and leads to the matrix equation below

M(r)TM(r)+βR(r)TR(r)α(r)=βR(r)TR(r)αB(r)+M(r)Tx(r)
(23)

The fact that the quantities in the above equation are all the functions of position vector r implies that they needs to be solved for each FEM node, as indicated before. Mathematically, Eq. (22) or Eq. (23) seeks the filtered least-squares solution to Eq. (17) on the basis of its smoothness around the baseline (the known background properties). The smoothing weight is controlled by the regularizing factor β, whose value is task-dependent. Again, the median filtering is applied to α(r) images as well as the resultant ημaf(r) and τ(r) ones, respectively, for a further suppression of the artifacts.

Fig. 1. FEM mesh and optode configuration employed in the study.

Table 1. Optical, fluorescent and geometrical parameters of Phantom 1

table-icon
View This Table
| View All Tables

Table 2. Optical, fluorescent and geometrical parameters of Phantom 2

table-icon
View This Table
| View All Tables

3. Validations

Fig. 2. (a) Original (top) and reconstructed (bottom) images of fluorescent yield and lifetime for Phantom 1, and (b) their profiles along X- and Y-axes.
Fig. 3. (a) Original (top) and reconstructed (bottom) images of fluorescent yield and lifetime for Phantom 2, and (b) original (red) and reconstructed (green) profiles along X-axis.
Fig. 4. (a) Schematic of the phantom used for the evaluating spatial resolution of the algorithm, (b) reconstructed images for CCS=13 mm (top), 15 mm (middle) and 17 mm (bottom), respectively, and (c) profiles of the reconstructed yield and lifetime images along the X-axis.

3.1 One-component case

In this case, we uniquely employ a pair of real transform-factors: p 1,2 = ∓0.1P, where P = 1/[1|μ(axB) c + 1/μam(B) c + τ(B)], for the Laplace transforms of the time-domain forward and inverse models, and a circulation number M = 20 as well as a relaxation parameter λ = 0.5 for the ART solution to Eq. (8), the derived linear system. The validation is firstly enforced by two numerical phantoms, for which the noiseless sets of the data-type Γ(p) are generated from the forward model with the same mesh as that to be used in the inverse model, to evaluate the intrinsic performance of the algorithm. In Phantom 1, four fluorescent circular targets with the same radius but different fluorescent properties are included to investigate the ability of the algorithm to discern the difference in the fluorescent properties of the targets, while Phantom 2 embeds two disks with same fluorescence parameters but different radiuses in order to probe the performance of the algorithm to reconstruct the target size. The optical properties of the targets in the two phantoms are the same as those of the background, which are set to μax,m = 0.035mm-1 and μ′sx,m = 1.0mm-1 for both the excitation and emission wavelengths. These values are in the range of the optical properties for in vivo muscle [31

31. F. Bevilacqua, D. Piguet, P. Marquet, J. D. Gross, B. J. Thromberg, and C. Depeursinge, “In vivo local determination of tissue optical properties: applications to human brain,” Appl. Opt. 38, 4939–4950 (1999). [CrossRef]

]. All the optical and fluorescent properties for Phantom 1 and 2 are listed in Table 1 and 2, respectively. The original and reconstructed images and their profiles along the axes for both the phantoms are shown in Fig. 2 and 3, respectively, where all the differences in the yield, lifetime and size among the targets are reasonably disclosed and further improvements in the spatial resolution and quantitative accuracy are desired. It is worthy to note in Fig. 3 that the reconstruction of the large sized target is much more quantitatively faithful than that of the small-sized one.

Secondly, we investigate the spatial resolution of the algorithm using a phantom with the same domain and background optical/fluorescent properties as the above examples. Two fluorescent disks with the same radius of 3 mm and optical/fluorescent properties of μax,m =0.035 mm-1 , μ′sx,m =1.0 mm-1, ημaf =0.005 mm-1 and τ = 1000 ps , but a variable center-to-center separation (CCS), are placed along the X-axis, as shown in Fig. 4(a). Figure 4(b) presents the reconstructed images and Fig. 4(c) their profiles along the X-axis, for three different values of the CCS: 13 mm, 15mm and 17 mm. It is seen from the results that the two disks can still be resolved as the CCS is less than 13 mm (the valleys drops from the peaks by about 9% and 5% for the yield and lifetime, respectively), and the recovery of the fluorescent yield appears to have better resolution but poorer quantitativeness (only about 45% of the true value is reconstructed by the peaks) than the lifetime in terms of the profiles.

Fig. 5. Investigation on the noise robustness of the algorithm by imaging the same phantom as in Fig. 4, with the CCS of the two target disks equals to 17 mm, for a varying SNR of 35 db (Top), 40 dB (Middle) and 45 dB (Bottom). (a) Reconstructed yield and lifetime images, and (b) their profiles along the X-axis.

3.2 Two-component case

Simultaneous reconstruction of two sets of fluorescent yield and lifetime images in the two-component case is in principle a direct extension to the one-component case but involves more complex mathematics, as aforementioned. Here we only present preliminary results for imaging a simple two-component phantom. The phantom has the same domain and background/target optical properties as in Fig. 4 and embeds two circular targets with the same radius of 4 mm and a fixed CCS of 25 mm. The background and targets all contain two kinds of fluorescent probes, as listed in Table 3. A set of 11 real transform-factors from -0.25P to +0.25P with a step of 0.05P is used for the Laplace transforms. This makes Eq. (17) overdetermined and ill-posed, with the condition number of its normal equation being 105 in the order of magnitude, for which the least-square solution can be approximately found with the Levenberg-Marquardt scheme. The reconstruction is performed with Gaussian noisy data of χ = 45 dB and a fixed regularizing factor β = 10-3. Figure 6 shows the original and reconstructed images for the two-component phantom, as well as their profiles along X-axis. It is observed that the reconstruction of the fluorescent yield correctly reflects the difference between the two components but both are underestimated, while the resultant images of the lifetime are somewhat distorted with the first component underestimated and the second overestimated.

Table 3. Fluorescent parameters in the two-component phantom

table-icon
View This Table
| View All Tables

4. Discussions and conclusions

In general sense, the Laplace transform of the time-resolved data in the above GPST-FMT algorithm can be implemented in the complex-domain. In our simulations only real transform-factors are employed. Compared to the methods using the Fourier transform, which is the imaginary-domain case of the Laplace transform and can be directly provided by a frequency-domain system, those that are based on the real-domain Laplace transform are less computationally intensive since no complex-number arithmetic is involved in the inversion procedure. Unlike the Fourier transform whose intensity and phase concern the information on the fluorescent yield and lifetime, respectively, the Laplace transform utilizes the information embedded in the full time-resolved data in an implicit way: the Laplace transform using a positive or negative transform-factor can be physically interpreted as an exponential time weighting that favors “early-time” or “late-time” features of the time-resolved curves that are jointly influenced by the fluorescent yield and the lifetime, and thus a pair of positive and/or negative frequencies would effectively distinguish between them. Originated from this idea, for the multi-component case with Nc targets, a proper choice of at least Nc transform-factor pairs would selectively pick up the signal weights that favor the different fluorescent emissions, and lead to a reliable separation among the fluorophores. A major disadvantages of the real-domain Lapalace transform is that it is now unclear if the information embedded in the time-resolved data can be extracted effectively, as with the Fourier transform, which mathematically complete along the imaginary domain.

Fig. 6. (a) Original and (b) reconstructed images of the two-component phantom with the fluorescence parameters listed in Table 3: The top row is for the first component and the bottom row for the second component in each sub-figure. The original (red) and reconstructed (green) profiles along X-axis are shown in (c) and (d) for Component 1 and Component 2 respectively.

The foundation of the above linear scheme of time-domain FMT lies on the assumption that the influence of the fluorophore absorption on the optical properties at both the excitation and emission wavelength is negligible, implicating that the targets to be imaged are of either low-contrast or small-size so that the evolvement of their fluorescent parameters during the reconstruction will not result in an evident change in the weighting matrix of Eq. (8). For accurately recovering targets of large-size and high-contrast, or for precisely quantifying fluorophores of low-concentration that is usual case in the diagnosis and pathogenesis research of early tumor, the influence of the probe absorption on the weighting matrix should be necessarily taken into account, even though some authors have demonstrated that the FMT image reconstruction might be insensitive to the background heterogeneity to some extend [32

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

]. This eventually leads to a nonlinear inverse issue that needs to be solved with an iterative scheme.

It is rather difficult to characterize the noise in the measurement since it can come from many sources. Generally, two types of noise sources are dominant in the measured signals. The first is due to thermal noise 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 Poisson distributed in the fluence rate. This model can be approximated with independent, but not identically distributed Gaussian random variables, with an expectation that fluence rate will need to be reasonably high to exceed the thermal noise floor of the detectors. Since the real Lapalce transform of the original time-resolved signal can be in broad sense regarded as a windowed intensity, we thus, for the most simplification, use the approximated shot noise model aforementioned. Because the noise model does not have a constant variance across the source-detector pairs, a whitening procedure might be incorporated for more robust reconstruction. This is accomplished by weighting Eq. (8) with the inverse covariance matirx

C1Γ(p)=C1W(p)x(p)
(24)

where C -1 = diag(σ(ξds,p)) is the inverse of the covariance matrix of the data-type. Eq. (24) requires, however, that the covariance matrix is either known as a prior knowledge or estimated as auxiliary parameters within an inversion procedure [5

05. A.B. Milstein, S. Oh, K.J. Webb, C.A Bouman, Q Zhang, D.A. Boas, and R.P. Millane, “Fluorescence optical diffusion tomography,” Appl. Opt. 42, 3081–94 (2003). [CrossRef] [PubMed]

]. It should be noted that the assumption of the constant SNR throughout the detectors is somewhat inconsistent with the real scenario where the SNR varies proportionally with the square root of the fluence rate present at each detector [5

05. A.B. Milstein, S. Oh, K.J. Webb, C.A Bouman, Q Zhang, D.A. Boas, and R.P. Millane, “Fluorescence optical diffusion tomography,” Appl. Opt. 42, 3081–94 (2003). [CrossRef] [PubMed]

]. Nevertheless, the shot-noise model with a constant SNR (most favorably set to the minimum) helps reasonably identify the lower limit of the noise-robustness of the algorithm since a time-resolved multi-channel system usually try to equalize the measured intensities through attenuator adjustment. According to the status of the time-resolved instrumentation, in particular a time-correlated single photon counting technique, a noise-robustness of SNR around 40 dB would be favored in practice.

The above examples of the image reconstruction are based on the absolute values of the Laplace-transformed time-varying flux. The evaluations show that the noise-robustness of this absolute imaging scheme is moderate. This indicates some potential difficulties in achieving a high-sensitive FMT for small animal, where SNR might be rather low due to the strong background fluorescence, the inevitable systematic noises as well as probably inaccurate estimation of the background optical and fluorescent properties. Fortunately, the differential imaging scheme is available in the most practice of FMT, where the measurements before and after intravascular injection of the fluorescent probes can be obtained and the difference between them is used to form the data-types that is to be used for the reconstruction, i.e.,

Γm(ξd,ζs,p)=[Γm(A)(ξd,ζs,p)Γm(B)(ξd,ζs,p)]+Γm(M)(ξd,ζs,p)
(25)

where Γm(B)d, ζs, p) and Γm(A)d, ζs, p) are the data-types before and after probe injection, respectively; Γm(M)d, ζs, p) the model prediction in terms of the background parameters. Differential imaging has been proved to be able to cancel many systematic and measuring noises and greatly improve the image quality [23–27

23. E.M.C. Hillman, J.C. Hebden, M. Schweiger, H. Dehghani, F.E.W. Schmidt, D.T. Delpy, and S.R. Arridge, “Time resolved optical tomography of the human forearm,” Phys. Med. Biol. 46, 1117–1130 (2002). [CrossRef]

]. It is also important to note that with the differential imaging scheme the three-dimensional FMT task can be performed in a layered imaging mode where the source-detector pairs are deployed to obtain the image at each of the multiple planes, by use of a two-dimensional inverse model. The validity of this method has been justified in DOT and FMT [10

10. V. Ntziachristos, C-H Tung, C. Bremer, and R. Weissleder, “Fluorescence molecular tomography resolves protease activity in vivo,” Nat. Med. 8, 757–60 (2002). [CrossRef] [PubMed]

, 26

26. Huijuan Zhao, Gao Feng, Yukari Tanikawa, Kazuhiro Homma, and Yukio 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]

].

In the two-component simulation, the regularized Levenberg-Marquardt method, i.e., Eq. (23), has been promisingly applied for generating the reliable α -images from the overdetermined but ill-conditioned linear system that are derived from the x -images at a number of different tramsform-factors, with a constant regularizing-factor β. Furthermore, it has been shown that an L-curve criterion for determining β can better trade off between the data misfit ||x(r) - M(r)α(r)|| and the solution-to-baseline seminorm ||R(r)[α(r) - αB(r)]||, but at a considerably increasing cost of computation. By comparison, we found that this regularized overestimation scheme can produce much higher quality of images than the simple matrix inversion using only four transform-factors. It is also pointed out here that, although we confine the fluorescent heterogeneities for the two components to the same regions for simplicity, it is unnecessary in the methodology for the fluorescence regions of the multiple probes to be the same as each other. The allowance for the maximal freedom of the probe distributions agrees with the fact that the multiple probes specifically bound to their respective biochemical molecules have little probability of the same fluorescent emission behavior, due to the intrinsic differences among the distributions of the molecule concentrations and among the efficacies of the fluorescent probes.

The Laplace transformed time-dependent diffusion model with a real transform-factor p, i.e., Eq. (1), can mathematically interpreted as the time-independent coupled diffusion equations with absorption coefficients of μav + p/c (v ∈ {x, m}). Therefore, together with considering the exponentially decaying response of the fluorescent emission, the choice of the transform-factor p should meet the condition that p ≥ max{- μav c, - 1/τ} so that the time-independent diffusion model is physically meaningful as well mathematically stable. Besides this criterion, there is no theoretical guide to the choice of the transform-factor pairs. Therefore the transform-factor pairs used in the above examples are somewhat empirical. Nevertheless, the above numerical validations of the methodology with these empirically chosen transform-factor pairs have achieved our goal with considerable success. The selection of the transform-factor pairs that can optimally characterize the original time-resolved signal will be a challenge in the future investigation.

With the objective of developing a general methodology for the image reconstruction of TD-FMT that can work with a broad range of fluorescent properties of the probes that might be encountered in applications, the fluorescence parameters in the above examples do not all exactly conforms to those of the common fluorescent dyes, such as Indocyanine green (ICG) (with the extinction coefficient of about 0.013 mm-1μM-1, the lifetime of about 0.56 ns and the quantum efficiency of 0.016 at the peak excitation wavelength of 778 nm) and Cy5.5 (with the extinction coefficient of about 0.019 mm-1μM-1, the lifetime of about 1 ns and the quantum efficiency of 0.23 at the peak excitation wavelength of 670 nm). Nevertheless, with regard to these two dyes, the fluorophore concentrations corresponding to the chosen fluorescent yields of the targets are all on the orders of 1μM and 100 nM, respectively, as preferably expected in vivo.

According to the above examples, the target size and contrast significantly influences the image quality. A formal investigation on these effects is shown in Fig. 7, where an increasing improvement in the quantitativeness is observed with the size increasing and/or the contrast decreasing. This observation conforms to that found in DOT. For the two-component case, because the lifetime values of the two components are reciprocal in Eq. (16), an underestimation of one lifetime value may result in an overestimation of another one, and vice versa. This adverse effect has been clearly seen in Fig. 5, where the lifetime reconstruction for the first component is underestimated while that for the second slightly overestimated. In this sense, an effective enhancement of the solution accuracy of the first-stage equation, i.e., Eq. (17), is critical to the reconstruction fidelity. Normally, this would be attainable through regularizing ill-posedness of the inverse issue and incorporating a priori knowledge. It is also argued that two-component reconstruction is target-dependent and would be improved for targets of large-size and low-contrast.

Fig. 7. Effects of target size (a) and contrast (b) on the reconstruction quantitativeness. (a) Target contrast fixed at 3:1 for a varying radius; (b) Target radius fixed at 4 mm for varying contrast. The same phantom as in Fig. 4 is used with the target CCS=25 mm, and the results are shown for the peak values in the reconstructed images.

It is clearly seen from the profiles of the reconstructed images, that the centers of the reconstructed targets appear to be deeper and the reconstructed sizes larger than the true ones. The latter is the intrinsic nature of diffuse light imaging modality that stems from ill-posedness of the inverse issue; while the reason for the former observation is complex and can be preliminarily ascribed to the use of the row-fashioned ART inversion. It is expected that both the adversities would be substantially suppressed with adoption of the whole-weight-matrix based inversion scheme and incorporation of the prior-knowledge based constrains.

In summary, we have presented a linear GPST methodology of TD-FMT for simultaneous reconstruction of fluorescent yields and lifetimes of multi-components, from which the algorithms for one- and two-component recovery are exemplified, respectively. The simulative validations of the algorithm for a 2-D domain have been performed for its abilities to discern differences in target size and fluorescence parameters, as well as for its noise-robustness and spatial resolution. The results have proved the effectiveness of the methodology. Also, the crucial issues on the implementation and validation of the algorithm have been discussed. We will explore phantom and in vivo experimental validations of the methodology in the on-going work.

References and links

01.

D.Y. Patthankar, A.U. Chen, B.W. Pogue, M.S. Patterson, and E.M. Sevick-Muraca, “Imaging of fluorescent yield and lifetime from multiply scattered light reemitted from random media,” Appl. Opt. 36, 2260–2272 (1997). [CrossRef]

02.

H. Jiang, “Frequency-domain fluorescent diffusion tomography: a finite-element-based algorithm and simulation,” Appl. Opt. 37, 5337–5343 (1998). [CrossRef]

03.

E.M. Sevick-Muraca, J.P. Houston, and M. Gurfinkel, “Fluorescence-enhanced, near infrared diagnostic imaging with contrast agent,” Curr. Opin. Chem. Biol. 6, 642–50 (2002). [CrossRef] [PubMed]

04.

E.J. Eppstein, D.J. Hawrysz, A. Godavarty, and E.M. Sevick-Muraca, “Three-dimensional, Bayesian image reconstruction from sparse and noisy data sets: Near-infrared fluorescence tomography,” Proc. Acad. Sci. Am. 99, 9619–9624 (2002). [CrossRef]

05.

A.B. Milstein, S. Oh, K.J. Webb, C.A Bouman, Q Zhang, D.A. Boas, and R.P. Millane, “Fluorescence optical diffusion tomography,” Appl. Opt. 42, 3081–94 (2003). [CrossRef] [PubMed]

06.

K. Licha, “Contrast agents for optical imaging,” Topics in Current Chemistry 222, 1–29 (2002). [CrossRef]

07.

Achilefu, R. Dorshow, J. Bugaj, and R. Rajapopalan, “Novel receptor-targeted fluorescent contrast agents for in vivo tumor imaging,” Invest. Radiol. 35, 479–485 (2000). [CrossRef] [PubMed]

08.

R. Weissleder, C.H. Tung, U. Mahmood, and A. Bogdanov, “In vivo imaging with protease-activated near-infrared fluorescent probes,” Nat. Biotechnol. 17, 375–378 (1999). [CrossRef] [PubMed]

09.

R.E. Campbell, O. Tour, A.E. Palmer, P.A. Steinbach, G.S. Baird, D.A. Zacharias, and R. Tsien, “A monometric red fluorescent protein,” Proc. Natl. Acad. Sci. USA 99, 7877–7882 (2002). [CrossRef] [PubMed]

10.

V. Ntziachristos, C-H Tung, C. Bremer, and R. Weissleder, “Fluorescence molecular tomography resolves protease activity in vivo,” Nat. Med. 8, 757–60 (2002). [CrossRef] [PubMed]

11.

V. Ntziachristos, C. Bremer, E.E. Graves, J. Ripoll, and R. Weissleder, “In vivo tomographic imaging of near-infrared fluorescent probes,” Molecular Imaging 1, 82–88 (2002). [CrossRef]

12.

S. Lam, F. Lesage, and X. Intes, “Time domain fluorescent diffuse optical tomography: analytical expressions,” Opt. Express 13, 2263–2275 (2005). [CrossRef] [PubMed]

13.

A.T.N. Kumar, J. Skoch, B.J. Bacskai, D.A. Boas, and A.K. Dunn, “Fluorescent-lifetime-based tomography for turbid media,” Opt. Lett. 30, 3347–3349 (2005). [CrossRef]

14.

X. Cong and G. Wang, “A finite-element-based reconstruction method for 3D fluorescence tomography,” Opt. Express 13, 9847–9857 (2005). [CrossRef] [PubMed]

15.

S.R. Cherry, “In vivo molecular and genomic imaging: new challenges for imaging physics,” Phys. Med. Biol. 49, R13–48 (2004). [CrossRef] [PubMed]

16.

T.F. Massoud and S.S. Gambhir, “Molecular imaging in living subjects: seeing fundamental biological processes in a new light,” Genes Dev. 17, 545–580 (2003). [CrossRef] [PubMed]

17.

A.D. Klose, V. Ntziahristos, and A.H. Hielschler, “The inverse source problem based on the reative trabsfer equation in optical molecular imaging,” J. Comput. Phys. 202, 323–345 (2002). [CrossRef]

18.

S.R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15, R41–93 (1999). [CrossRef]

19.

F. Gao, H. Zhao, and Y. Yamada, “Improvement of image quality in diffuse optical tomography by use of full time-resolved data,” Appl. Opt. 41, 778–791 (2002). [CrossRef] [PubMed]

20.

R. Model, M. Orlt, and M. Walzel, “Reconstruction algorithm for near-infrared imaging in turbid media by means of time-domain data,” J. Opt. Soc. Am. A 14, 313–323 (1997). [CrossRef]

21.

M. Schweiger and S.R. Arridge, “Application of temporal filters to time resolved data in optical tomography,” Phys. Med. Biol. 44, 1699–1717 (1999). [CrossRef] [PubMed]

22.

F. Gao, P. Poulet, and Y. Yamada, “Simultaneous mapping of absorption and scattering coefficients from a three-dimensional model of time-resolved optical tomography,” App. Opt 39, 5898–5910 (2001). [CrossRef]

23.

E.M.C. Hillman, J.C. Hebden, M. Schweiger, H. Dehghani, F.E.W. Schmidt, D.T. Delpy, and S.R. Arridge, “Time resolved optical tomography of the human forearm,” Phys. Med. Biol. 46, 1117–1130 (2002). [CrossRef]

24.

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]

25.

F. Gao, H. Zhao, Y. Tanikawa, and Y. Yamada, “Optical tomographic mapping of cerebral haemodynamics by time-domain detection: methodology and phantom validation,” Phys. Med. Biol. 49, 1055–1078 (2004). [CrossRef] [PubMed]

26.

Huijuan Zhao, Gao Feng, Yukari Tanikawa, Kazuhiro Homma, and Yukio 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]

27.

F. Gao, Y. Tanikawa, H.J. Zhao, and Y. Yamada, “Semi-three-dimensional algorithm for time-resolved diffuse optical tomography by use of the generalized pulse spectrum technique,” Appl. Opt. 41, 7346–7358 (2002). [CrossRef] [PubMed]

28.

W.G. Egan and T.W. Hilgeman, Optical Properties of Inhomogeneous Materials, (Academic, New York, 1979).

29.

A.C. Kak and M. Slaney, Principle of Computerized Tomographic Imaging, (IEEE Press, New York, 1988).

30.

F. Gao, H. Niu, H. Zhao, and H. Zhang, “The forward and inverse models in time-resolved optical tomography imaging and their finite-element method solutions,” Image and Vision Computing 16, 703–712 (1998). [CrossRef]

31.

F. Bevilacqua, D. Piguet, P. Marquet, J. D. Gross, B. J. Thromberg, and C. Depeursinge, “In vivo local determination of tissue optical properties: applications to human brain,” Appl. Opt. 38, 4939–4950 (1999). [CrossRef]

32.

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]

OCIS Codes
(170.3010) Medical optics and biotechnology : Image reconstruction techniques
(170.3660) Medical optics and biotechnology : Light propagation in tissues
(170.3880) Medical optics and biotechnology : Medical and biological imaging
(170.6280) Medical optics and biotechnology : Spectroscopy, fluorescence and luminescence
(170.6920) Medical optics and biotechnology : Time-resolved imaging
(170.6960) Medical optics and biotechnology : Tomography

ToC Category:
Medical Optics and Biotechnology

History
Original Manuscript: May 17, 2006
Revised Manuscript: July 24, 2006
Manuscript Accepted: July 24, 2006
Published: August 7, 2006

Virtual Issues
Vol. 1, Iss. 9 Virtual Journal for Biomedical Optics

Citation
Feng Gao, Huijuan Zhao, Yukari Tanikawa, and Yukio 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/oe/abstract.cfm?URI=oe-14-16-7109


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. D.Y. Patthankar, A.U. Chen, B.W. Pogue, M.S. Patterson, and E.M. Sevick-Muraca, "Imaging of fluorescent yield and lifetime from multiply scattered light reemitted from random media," Appl. Opt. 36, 2260-2272 (1997). [CrossRef]
  2. H. Jiang, "Frequency-domain fluorescent diffusion tomography: a finite-element-based algorithm and simulation," Appl. Opt. 37, 5337-5343 (1998). [CrossRef]
  3. E.M. Sevick-Muraca, J.P. Houston, and M. Gurfinkel, "Fluorescence-enhanced, near infrared diagnostic imaging with contrast agent," Curr. Opin. Chem. Biol. 6, 642-50 (2002). [CrossRef] [PubMed]
  4. E.J. Eppstein, D.J. Hawrysz, A. Godavarty and E.M. Sevick-Muraca, "Three-dimensional, Bayesian image reconstruction from sparse and noisy data sets: Near-infrared fluorescence tomography," Proc. Acad. Sci. Am. 99, 9619-9624 (2002). [CrossRef]
  5. A.B. Milstein, S. Oh, K.J. Webb, C.A. Bouman, Q Zhang, D.A. Boas, and R.P. Millane, "Fluorescence optical diffusion tomography," Appl. Opt. 42, 3081-94 (2003). [CrossRef] [PubMed]
  6. K. Licha, "Contrast agents for optical imaging," Topics in Current Chemistry 222, 1-29 (2002). [CrossRef]
  7. Achilefu, R. Dorshow, J. Bugaj and R. Rajapopalan, "Novel receptor-targeted fluorescent contrast agents for in vivo tumor imaging," Invest. Radiol. 35, 479-485 (2000). [CrossRef] [PubMed]
  8. R. Weissleder, C.H. Tung, U. Mahmood and A. Bogdanov, "In vivo imaging with protease-activated near-infrared fluorescent probes," Nat. Biotechnol. 17, 375-378 (1999). [CrossRef] [PubMed]
  9. R.E. Campbell, O. Tour, A.E. Palmer, P.A. Steinbach, G.S. Baird, D.A. Zacharias and R. Tsien, "A monometric red fluorescent protein," Proc. Natl. Acad. Sci. USA 99, 7877-7882 (2002). [CrossRef] [PubMed]
  10. V. Ntziachristos, C-H Tung, C. Bremer and R. Weissleder, "Fluorescence molecular tomography resolves protease activity in vivo," Nat. Med. 8, 757-60 (2002). [CrossRef] [PubMed]
  11. V. Ntziachristos, C. Bremer, E.E. Graves, J. Ripoll, and R. Weissleder, "In vivo tomographic imaging of near-infrared fluorescent probes," Molecular Imaging 1, 82-88 (2002). [CrossRef]
  12. S. Lam, F. Lesage, and X. Intes, "Time domain fluorescent diffuse optical tomography: analytical expressions," Opt. Express 13, 2263-2275 (2005). [CrossRef] [PubMed]
  13. A.T.N. Kumar, J. Skoch, B.J. Bacskai, D.A. Boas, and A.K. Dunn, "Fluorescent-lifetime-based tomography for turbid media," Opt. Lett. 30, 3347-3349 (2005). [CrossRef]
  14. X. Cong and G. Wang, "A finite-element-based reconstruction method for 3D fluorescence tomography," Opt. Express 13, 9847-9857 (2005). [CrossRef] [PubMed]
  15. S.R. Cherry, "In vivo molecular and genomic imaging: new challenges for imaging physics," Phys. Med. Biol. 49, R13-48 (2004). [CrossRef] [PubMed]
  16. T.F. Massoud and S.S. Gambhir, "Molecular imaging in living subjects: seeing fundamental biological processes in a new light," Genes Dev. 17, 545-580 (2003). [CrossRef] [PubMed]
  17. A.D. Klose, V. Ntziahristos, and A.H. Hielschler, "The inverse source problem based on the reative trabsfer equation in optical molecular imaging," J. Comput. Phys. 202, 323-345 (2002). [CrossRef]
  18. S.R. Arridge, "Optical tomography in medical imaging," Inverse Probl. 15, R41-93 (1999). [CrossRef]
  19. F. Gao, H. Zhao, and Y. Yamada, "Improvement of image quality in diffuse optical tomography by use of full time-resolved data," Appl. Opt. 41, 778-791 (2002). [CrossRef] [PubMed]
  20. R. Model, M. Orlt, and M. Walzel, "Reconstruction algorithm for near-infrared imaging in turbid media by means of time-domain data," J. Opt. Soc. Am. A 14, 313-323 (1997). [CrossRef]
  21. M. Schweiger and S.R. Arridge, "Application of temporal filters to time resolved data in optical tomography," Phys. Med. Biol. 44, 1699-1717 (1999). [CrossRef] [PubMed]
  22. F. Gao, P. Poulet and Y. Yamada, "Simultaneous mapping of absorption and scattering coefficients from a three-dimensional model of time-resolved optical tomography," App.Opt 39, 5898-5910 (2001). [CrossRef]
  23. E.M.C. Hillman, J.C. Hebden, M. Schweiger, H. Dehghani, F.E.W. Schmidt, D.T. Delpy and S.R. Arridge, "Time resolved optical tomography of the human forearm," Phys. Med. Biol. 46, 1117-1130 (2002). [CrossRef]
  24. 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]
  25. F. Gao, H. Zhao, Y. Tanikawa, and Y. Yamada, "Optical tomographic mapping of cerebral haemodynamics by time-domain detection: methodology and phantom validation," Phys. Med. Biol. 49, 1055-1078 (2004). [CrossRef] [PubMed]
  26. Huijuan Zhao, Feng Gao, Yukari Tanikawa, Kazuhiro Homma, and Yukio 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]
  27. F. Gao, Y. Tanikawa, H.J. Zhao and Y. Yamada, "Semi-three-dimensional algorithm for time-resolved diffuse optical tomography by use of the generalized pulse spectrum technique," Appl. Opt. 41, 7346-7358 (2002). [CrossRef] [PubMed]
  28. W.G. Egan and T.W. Hilgeman, Optical Properties of Inhomogeneous Materials, (Academic, New York, 1979).
  29. A.C. Kak and M. Slaney, Principle of Computerized Tomographic Imaging, (IEEE Press, New York, 1988).
  30. F. Gao, H. Niu, H. Zhao and H. Zhang, "The forward and inverse models in time-resolved optical tomography imaging and their finite-element method solutions," Image and Vision Computing 16, 703-712 (1998). [CrossRef]
  31. F. Bevilacqua, D. Piguet, P. Marquet, J. D. Gross, B. J. Thromberg and C. Depeursinge, "In vivo local determination of tissue optical properties: applications to human brain," Appl. Opt. 38, 4939-4950 (1999). [CrossRef]
  32. 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]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.


« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited