Open Access
Add to BibSonomyAbstract
We present a new method allowing the reconstruction of 3D time-domain diffuse optical tomography images, based on the time-dependent diffusion equation and an iterative algorithm (ART) using specific points on the temporal profiles. The first advantage of our method versus the full time-resolved scheme consists in considerably reducing the inverse problem resolution time. Secondly, in the presence of scattering heterogeneities, our method provides images of better quality comparatively to classical methods using full-time data or the first moments of the profiles.
©2011 Optical Society of America
1. Introduction
Diffuse optical tomography (DOT) based on the Near-Infrared (NIR) spectral window is a rapidly growing field with applications ranging from human breast cancer detection [1] and brain functional imaging [2], to small animals imaging [3]. In the therapeutic window [4] where absorption and scattering are low, it is possible to detect photons that have traveled several centimeters through biological tissue [5]. Thus, coupled with an accurate photon migration model, DOT enables mapping of the optical properties of deep tissues from boundary measurements. The major interest of this technique is the use of low-cost instrumentation and the ability to monitor oxy- and deoxy-hemoglobin concentrations [6–8].
Over the past decade, significant progress has been achieved in three main DOT measurements techniques: 1) continuous wave (CW) systems, measuring the attenuation of CW light traveling through a given object; 2) frequency-domain (FD) systems using the phase shift and modulation depth between the irradiating and the detected signals; and 3) time-resolved (TR) systems measuring a time-dependent flux at the boundary of a studied object after its excitation by an ultrashort pulsed laser [9–12].
Alternatively, time-weighted moments derived from Mellin transform of the original TR measurements can be used to reduce the large amount of data [13–17]. This method considerably accelerates the image reconstruction process since it uses time-weighted moments predicted from a time-dependent diffusion equation (TDE) without explicitly resolving the TR profiles. These types of data are often preferred due to their robustness to noise and model mismatch [14]. However, this type of approach reduces the amount of information used in the reconstruction algorithm and can be considered as equivalent to FD-DOT when using the first 2 moments, energy and mean time of flight.
However DOT image reconstruction algorithms, using time-weighted moments, suffer from poor spatial resolution and quantitativeness [14–18], despite the suppression of crosstalk between the absorption and scattering spaces, which can be achieved by using a proper combination of time-weighted moment types [14,15]. Therefore, to overcome the lack of accuracy of such reconstruction algorithms, the exploitation of entire profiles or some of their points should be considered taking in mind the sensitivity to counting noise and time uncertainties [19–23].
In this paper, we present a new method based on the time-dependent diffusion equation and an iterative algorithm (ART) using specific points on the temporal profiles. Our method considerably accelerates the reconstruction process compared to the full TR scheme. Moreover, in the presence of scattering inhomogeneities, our method provides better quality images whereas the classical methods using full time or time-weighted moments data fail to discriminate between absorbing and scattering heterogeneities.
2. Method - the studied phantom
A numerical cylinder phantom was used, in order to test our algorithm [24]. Usually, the presence of scattering heterogeneity blurs the reconstructed absorption image. Although absorption is a major concern in physiological studies, the simultaneous mapping of both parameters is required to obtain the genuine absorption distribution. The model geometry was as follows: C, cylinder 60 mm high, 28 mm in diameter, centered on x = y = z = 0. The optical background properties of the model were: n = 1.54, µa = 0.005 mm−1, µs’ = 0.6 mm−1. The cylinder contained 2 rods, 5 mm in diameter and centered on:
- C1 at position x = 0 mm, y = 0 mm, z = 7 mm,
- C2 at position x = 0 mm, y = 0 mm, z = −7 mm,
Investigations were performed on two objects having the same geometry described above with different optical properties. First, object_1 contained only absorbing heterogeneities. Cylinders C1 and C2 were chosen to be 1.5 and twice as absorbing as the background, respectively. Second, object_2 contained scattering and absorbing heterogeneities, cylinder C1 was chosen to be 1.25 times as scattering and C2 twice as absorbing as the background. Figure 1 represents the geometry and mesh of the phantom described above. The mesh shown on Fig. 1 contained 3812 nodes, 18500 tetrahedral elements.
3. DOT forward model
3.1 Photon transport model
Image reconstruction techniques, in DOT, require a prevalent description of photon migration through the studied object. Using time-resolved data, this step consists in computing the temporal profiles Γ, or TPSF (Time Point Spread Function), of photons coming out of the surface of an object with known optical properties, for a given pulsed laser excitation.
It is generally accepted that with some assumptions the diffusion equation, i.e., the P1 approximation to the radiative transfer equation, is an efficient model of light transport in scattering media [20]. The time-dependent DE approximation expressed in terms of the photon density u(r,t) is given by [25–27]
where u(r,t) is the photon density at position r and time t, µa, the absorption coefficient and κ the diffusion coefficient, κ = c/3(μa + μs’) with µs’ being the reduced scattering coefficient and c the speed of light in the medium. δ(t, r0) is the Dirac pulse modeling the isotropic sources positioned under the irradiated surface of the studied object, at a distance from the surface imposed by the approximation to DE and equal to 1/µs’.We used the Neumann boundary condition for the diffusion equation:
with q being a function of the reflection coefficient at the surface, and the vector normal to the surface of the object at position r [27].In this study, irradiation and measurements were performed with a non-contact set-up. In fact the extremities of the optical fibers, used as sources and detectors, were mounted on a cylinder (Ø = 40 mm) surrounding the object of smaller diameter (Ø = 28 mm). In the gap between the fibers and the surface of the object, photons propagate rectilinearly delineating a cone-like pattern, with the fibers extremities as apexes and angles depending on their numerical apertures [24].
3.2 Generation of simulated measurement data sets
Simulated measurement data sets were computed for the model in Eq. (2) with a given distribution of {μa,μ’s} and using a mesh finer than the one described in section 2 (23242 nodes, 149545 tetrahedral elements). The data set consisted of the theoretical temporal profiles Γs,d(t) obtained using the Comsol Multiphysics® time-dependent solver for the 112 source-detector pairs obtained with 16 sources and 7 detectors per source at the model surface.
Considering the stochastic noise associated with TCSPC-based MCP-PMT detection, a Poisson noise was added to the simulated data [28]. Then the stochastic noise at any point of the temporal profile was associated with the SNR as
with N(t) the number of photon counts at time t.The time-weighted moments were then obtained for each source-detector pair:
The total intensity:
The mean time of flight of photon:
4. DOT inverse problem
The inverse problem consists in mapping one or more internal optical characteristics from a series of boundary measurements. This, considering our TR-DOT configuration, could be rephrased as computing the distribution {μa,μs’} from the measured TPSF, for a given distribution of light sources and detectors.
The Algebraic Reconstruction Technique (ART) was used to achieve the best fit between the measured and the computed quantities by selecting a proper distribution of the absorption and the scattering coefficients [29]. The estimated values x = [μa,μ’s] were iteratively updated using:
with J(D) being the Jacobian matrix corresponding to the data type D, and λ the relaxation parameter proven in the range {0,2} [15,28]. At each iteration step, the Δx were computed on all the used data D and added. Methods using three different data types D are presented below:4.1 Method 1: optimization using the time-weighted moments
The inverse problem was first solved by finding the best fit between the time-weighted moments (total intensityand meantime of flight of photons) of the measured and simulated TPSF. The advantages of this approach, which consists in applying the Mellin transform allowing to predict time-weighted moments without explicit computation of the temporal profiles, have been largely described by Arridge et al [13,14].
4.2 Method 2: optimization using a large number of points of the TPSF
It has been largely investigated and demonstrated that using the entire set of points of temporal profiles in a fully time-resolved reconstruction algorithm provides the most accurate estimation of the optical properties of the studied medium but this method is very time-consuming [29]. Nevertheless, when performed on some points of the temporal profiles, the optimization process time can be substantially reduced. Figure 2 (left) shows the 19 points used to compute both the absorption and reduced scattering coefficients. The central point was set at the maximum profile amplitude and the other 18 were distributed half and half on the rising edge and tail of the profile between 10% and 90% of the maximum amplitude by 10% increments.
Fig. 2 Choice of the points used in the optimization process according to method 2 (left) and method 3 (right). Methods described in section 4.2 and 4.3.
4.3 Method 3: optimization using a reduced number of points of the TPSF
In this method, we present a new way to choose the points on the temporal profiles. Optimization is performed using different times, upon updating of the absorption coefficient (μa_new = μa_old + ∆μa) or the reduced scattering coefficient (μs’_new = μs’_old + ∆μs’). This method is inspired from the idea that the contrast obtained by late gating of the TPSF is more sensitive to absorption heterogeneities than early gating contrast [30–33] and was recently used by Venugopal et al. [23]. In order to justify this method, we compared three TPSFi (i = 1,2,3), simulated on a homogeneous cylinder (Ø = 28 mm) with the source positioned at x = y = 0 mm, z = 14 mm and detector positioned at x = y = 0 mm, z = −14 mm (180° from the source position).
Normalization of the maximum amplitudes of three temporal profiles TPSFi (i = 1, 2, 3) computed using different {µai,µsi’} with µa2 ≠ µa1 = µa3, and µs3’ ≠ µs1’ = µs2’, enable point-to-point comparisons. Figure 2 (left) shows the influence of the TPSF due to an increase of 50% of μa or μs’. Indeed, using TPSF1 as a reference, one can observe that varying μs’ causes the rising edge of TPSF3 to shift more than that of TPSF2. Information about a change in µa in the medium is better extracted from the tail of the TPSF, since presented with logarithm scale, the slope of this tail is independent from µs’ and proportional to µa only.
Accordingly, the points used for the optimization process were chosen on the rising edge of the TPSF or on its tail to compute Δμs’ or Δμa, respectively. This method reduces the number of considered points by 81% compared to method 2. In fact, the inverse problem is solved by using only 3 points chosen on the rising edge at 20%, 50% and 80% of the maximum amplitude of the profile when reconstructing µs’ images, and only 4 points on the profile tail, selected at 10%, 20%, 50% and 80% of the maximum amplitude, when reconstructing µa images, as shown in Fig. 2 (right).
The choice of these points is mainly heuristics. The use of points on the rising edge when reconstructing µs’, and points on the tail of the TPSF when reconstructing µa can also be justified with a theoretical model of light propagation. We used the Patterson equation R(t,µa,µs’) for an infinite medium [34] and demonstrated that its partial derivative with respect to µa, multiplied by µa, (∂R/∂µa).µa, is equal to its partial derivative with respect to µs’, multiplied by µs’, (∂R/∂µs’).µs’, at a time which is always very close from that of the maximum of the TPSF. The choice of the amplitude of the individual points was mainly done by considering their time separation and the signal to noise ratio to get valuable information. The optimal number of points was experimentally found to be seven. Using more points did not improve significantly the quality of the reconstructed images but using fewer points degrades it drastically.
5. Results and discussion
5.1 Simulated data sets
Simulated data sets were computed by solving Eq. (2) using the FEM Comsol solver with 5 ps time steps and the mesh described in subsection 3.2, finer than that used in the optimization step, presented in section 2. Figure 3 (left) shows the diffused photon-density map inside object_2 in the fibers plane (x = 0), computed 1 ns after illumination. Fiber F1 is used as the source. The light beam entered object_2 at r0 (x = z = 0 mm, y = −14 mm), corresponding to an isotropic source positioned at 1/µs’ under r0: (x = z = 0 mm and y = −12.33 mm).
Fig. 3 Photon-density map inside object_2 (slice at x = 0), 1ns after the light pulse and position of the 3 detection fibers at 135°, 180° and −135° from the exciting fiber (F1) (left), the temporal profiles detected (right).
The position of the detection fibers for which the computed TPSF are presented is illustrated in the same figure. The detection fiber at 135° from the source fiber F1 was positioned to face the rod C2 and the one at −135° was positioned to face the rod C1. As can be observed on the temporal profiles (Fig. 3 (right)) that the absorbing/scattering inhomogeneities are responsible for differences between the TPSFs detected at + 135° and at −135°. The simulated TPSFs were then normalised to those acquired using our TCSPC-based MCP-PMT detection system before adding the Poisson noise. We used typical count rates and acquisition times corresponding to a maximum of about 105 photons.
5.2 Reconstruction results
In order to reduce the computation time necessary to resolve the direct problem, the maximum time step was fixed at 20 ps, between zero and 3000 ps. The stop criterion was defined as follow: the algorithm was stopped if the 3 global errors of reconstruction computed at 3 successive iterations after the iteration (i) are bigger than the error computed at iteration i. Then, the best reconstruction was attributed to the iteration i. The maximum number of iterations was set to 25, if the stop criterion was not satisfied.
5.2.1 Phantom object_1
The results obtained using method 2, are supposed to be the most accurate since it uses the richest information data [22,27]. Actually, this rule is only true when the studied medium has a homogeneous scattering coefficient. Thus, Δµs’ is usually set to naught and Eq. (8) is used to compute Δµa only. However, in-vivo tissues exhibit heterogeneous scattering coefficients so that Δµs’ must be computed (Δµs’ ≠ 0).
Figure 4 shows the reconstruction result of the numerical phantom object_1 containing only absorbing heterogeneities, using method 2 with Δµs’ ≠ 0. We can clearly notice that the algorithm diverges. The reconstructed µa map is totally false and thus the position of the rods is not localized. The reconstructed µs’ presented on Fig. 4 (bottom right) shows that the algorithm is able to detect the position of two rods but is not able to estimate the correct values. However, the attenuation of light traveling through the object is not only interpreted as increased absorption but also as increased scattering (crosstalk) which distorts the absorption map (0,58 mm−1<µs'<0,76 mm−1). Hence, to obtain acceptable results, Δµs’ is set to naught when reconstructing objects which contain only absorbing heterogeneities.
Fig. 4 Column 1: true and reconstructed μa, column 2: true and reconstructed μs’, using method 2 with Δµs’ ≠ 0. Slice at x = 0.
Figure 5 shows the true and the reconstructed µa and µs’ maps obtained on the numerical phantom object_1 using the three reconstruction methods. According to the results presented on Fig. 5 (left), we notice that method 1 (0,58<µs’<0,61) and method 3 (0,59<µs’<0,6) are stable when reconstructing µs’ unlike method 2 (µs’ has to be null). Figure 5 (right) shows that the three methods yield almost the same result when reconstructing µa. As cited in subsection 4.3, method 3 allows to solve the inverse problem by using 7 points on the profiles. This method reduces reconstruction time by approximately 67% compared to method 2.
Fig. 5 Column 1: true and reconstructed μa, column 2: true and reconstructed μs’, using method 1, 2, 3 respectively. Slice at x = 0.
We used an explicit parameter Edata to compare reconstructed images and the true data in order to evaluate the three methods, defined by Eq. (9):
with P being the number of pixels in the region of interest (ROI). The ROI can be either background, target or bleed-through. In a given image type (µa or µs’) the target is defined as the region having that coefficient different from the background. The bleed-through is defined as the region with the other coefficient different from the background. The Edata parameter computed either on the background or the target is called error and the Edata computed on the bleed-through is called crosstalk.Concerning object_1, both cylinders C1 + C2 were considered as the targets and C-(C1 + C2) represented the background.
Table 1 lists the errors and crosstalk used to compare the reconstructed images using methods 1, 2 and 3.
Table 1. Errors and Crosstalk on the Reconstructed µa, µs’ Maps of Object_1, Slice at x=0
The three methods are almost equivalent when reconstructing the absorption image but they differ when considering the scattering image. Method 2 does not allow to reconstruct µs’ images and method 3 is more accurate than method 1. Figure 6 (left,right) shows the profiles along the z-axis of the reconstructed µa and μs’ presented on Fig. 5, respectively.
The profiles presented in Fig. 6 (left,right) show that when reconstructing object exhibiting a scattering homogeneity method 2 provides a slightly more precise estimation of µa than either method 1 or 3. The reconstructed µa using method 1 and 3 look similar but method 3 shows a better stability when reconstructing µs’. In fact, the small increase in the reconstructed µs’ is centered on the rods (C1,C2) owing to crosstalk.
5.2.2 Phantom object_2
Phantom object_2, containing one absorbing and one scattering rod, was used to model in-vivo tissues. In the following results, Δµs’ was not set to naught and computed using Eq. (8). The µa and µs’ maps reconstructed using method 1, 2 and 3 are presented in Fig. 7 .
Fig. 7 Column 1: true and reconstructed μa, column 2: true and reconstructed μs’, using method 1, 2, 3 respectively. Slice at x = 0.
The reconstructed absorption and scattering maps presented in Fig. 7 show that method 1 and 2 failed to separate the two targets (crosstalk). However, method 3 allowed a clear distinction between the two targets. Moreover, method 3 allowed to detect the scattering target and showed good background homogeneity equal to 0.6 mm−1.
In phantom object_2, the cylinder C1 is the target in the reconstructed µs’ images and the bleed-through in the µa images. The cylinder C2 is the target in the reconstructed µa images and the bleed-through in the µs’ images. The background is chosen to be C-(C1 + C2).
The errors and crosstalk used to compare the reconstructed images using methods 1, 2 and 3 are listed in Table 2 .
Table 2. Errors and Crosstalk on the Reconstructed µa, µs’ Maps of Object_2, Slice at x=0
Figure 8 (left,right) show the profiles along the z-axis of the reconstructed µa and μs’ maps, presented on Fig. 7, respectively. These profiles show clearly that methods 1 and 3 provided more accurate reconstructions than method 2, in both cases owing to the detection of the more scattering rod. In fact, at the first iteration, method 2 fostered absorption in both rods. Thus, light attenuation was interpreted as increased absorption rather than scattering and huge crosstalk between absorption and scattering was observed.
The profiles on Fig. 8 show that method 3 allowed the best distinction between the rods C1 and C2. Furthermore, the reconstructed background values were very stable in both reconstructions of µa and µs’.
6. Conclusion
We introduced a new method which, from only seven specifics points on the temporal profile, allows to reconstruct 3D time-domain diffuse optical tomography images, based on the time-dependent diffusion equation and an iterative algorithm (ART). The first advantage of our method versus the full time-resolved scheme consists in reducing the inverse problem resolution time by approximately 67%. Secondly, in presence of scattering heterogeneities, our method yields good quality images where classical methods, using points arbitrarily chosen on the temporal profile or its time-weighted moments, fail to do so owing to crosstalk between absorption and scattering.
This method can be especially interesting when reconstructing images from data acquired with an experimental setup based on a time-gated CCD camera [23,35]. The acquisition of a limited, optimized, number of data sets will reduce the measurement time and/or improve the time resolution.
These results obtained on the resolution of the inverse problem, and previous results on the resolution of the direct problem in experiments involving no-contact set-ups [24] should give us the opportunity to reconstruct 3D images of small animals using our laboratory-built TR-DOT equipment.
Both the experimental setup and forward modeling are ready to perform fluorescence imaging, but the inverse problem yet requires further work on the optimization process.
Acknowledgments
This work was supported by the Conseil Régional d'Alsace. The authors wish to thank Nathalie Heider and François Xavier Blé for thoroughly reading the manuscript.
References and links
1. B. J. Tromberg, B. W. Pogue, K. D. Paulsen, A. G. Yodh, D. A. Boas, and A. E. Cerussi, “Assessing the future of diffuse optical imaging technologies for breast cancer management,” Med. Phys. 35(6), 2443–2451 (2008). [CrossRef] [PubMed]
2. J. C. Hebden, “Advances in optical imaging of the newborn infant brain,” Psychophysiology 40(4), 501–510 (2003). [CrossRef] [PubMed]
3. V. Ntziachristos, J. Ripoll, L. V. Wang, and R. Weissleder, “Looking and listening to light: the evolution of whole-body photonic imaging,” Nat. Biotechnol. 23(3), 313–320 (2005). [CrossRef] [PubMed]
4. W. G. Egan and T. W. Hilgeman, Optical Properties of Inhomogeneous Materials (Academic, 1979).
5. A. G. Yodh and B. Chance, “Spectroscopy and imaging with diffusing light,” Phys. Today 48(3), 34–40 (1995). [CrossRef]
6. X. Intes and B. Chance, “Non-PET functional imaging techniques: optical,” Radiol. Clin. North Am. 43(1), 221–234 (2005). [CrossRef] [PubMed]
7. F. F. Jöbsis, “Noninvasive, infrared monitoring of cerebral and myocardial oxygen sufficiency and circulatory parameters,” Science 198(4323), 1264–1267 (1977). [CrossRef] [PubMed]
8. Y. Lin, G. Lech, S. Nioka, X. Intes, and B. Chance, “Noninvasive, low-noise, fast imaging of blood volume and deoxygenation changes in muscles using light-emitting diode continuous-wave imager,” Rev. Sci. Instrum. 73(8), 3065–3074 (2002). [CrossRef]
9. S. B. Colak, D. G. Papaioannou, G. W. ’t Hooft, M. B. van der Mark, H. Schomberg, J. C. Paasschens, J. B. Melissen, and N. A. van Asten, “Tomographic image reconstruction from optical projections in light-diffusing media,” Appl. Opt. 36(1), 180–213 (1997). [CrossRef] [PubMed]
10. H. Jiang, K. D. Paulsen, U. L. Osterberg, B. W. Pogue, and M. S. Patterson, “Optical image reconstruction using frequency-domain data: simulations and experiments,” J. Opt. Soc. Am. A 13(2), 253–266 (1996). [CrossRef]
11. M. Schweiger, S. R. Arridge, and D. T. Delpy, “Application of the finite-element method for the forward and inverse models in optical tomography,” J. Math. Imaging Vis. 3(3), 263–283 (1993). [CrossRef]
12. H. Eda, I. Oda, Y. Ito, Y. Wada, Y. Oikawa, Y. Tsunazawa, M. Takada, Y. Tsuchiya, Y. Yamashita, M. Oda, A. Sassaroli, Y. Yamada, and M. Tamura, “Multichannel time-resolved optical tomographic imaging system,” Rev. Sci. Instrum. 70(9), 3595–3602 (1999). [CrossRef]
13. S. R. Arridge, M. Schweiger, and D. T. Delpy, “Iterative reconstruction of near-infrared absorption images,” Proc. SPIE 1767, 372–383 (1992). [CrossRef]
14. M. Schweiger and S. R. Arridge, “Application of temporal filters to time resolved data in optical tomography,” Phys. Med. Biol. 44(7), 1699–1717 (1999). [CrossRef] [PubMed]
15. F. Gao, P. Poulet, and Y. Yamada, “Simultaneous mapping of absorption and scattering coefficients from a three-dimensional model of time-resolved optical tomography,” Appl. Opt. 39(31), 5898–5910 (2000). [CrossRef]
16. E. M. C. Hillman, J. C. Hebden, F. E. R. 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(9), 3415–3427 (2000). [CrossRef]
17. M. Torregrossa, C. V. Zint, and P. Poulet, “Effects of prior MRI information on image reconstruction in diffuse optical tomography,” Proc. SPIE 5143, 29–40 (2003). [CrossRef]
18. H. Dehghani, M. E. Eames, P. K. Yalavarthy, S. C. Davis, S. Srinivasan, C. M. Carpenter, B. W. Pogue, and K. D. Paulsen, “Near infrared optical tomography using NIRFAST: algorithm for numerical model and image reconstruction,” Commun. Numer. Methods Eng. 25(6), 711–732 (2009). [CrossRef]
19. R. L. Barbour, H. L. Graber, Y. Wang, J. H. Chang, and R. Aronson, “A perturbation approach for optical diffusion tomography using continuous-wave and time-resolved data,” Proc. SPIE IS11, 87–120 (1993).
20. A. H. Hielscher, A. D. Klose, and K. M. Hanson, “Gradient-based iterative image reconstruction scheme for time-resolved optical tomography,” IEEE Trans. Med. Imaging 18(3), 262–271 (1999). [CrossRef] [PubMed]
21. R. Model, M. Orlt, M. Walzel, and R. Hünlich, “Reconstruction algorithm for near-infrared imaging in turbid media by means of timedomain data,” J. Opt. Soc. Am. A 14(1), 313–323 (1997). [CrossRef]
22. 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(4), 778–791 (2002). [CrossRef] [PubMed]
23. V. Venugopal, J. Chen, and X. Intes, “Development of an optical imaging platform for functional imaging of small animals using widefield excitation,” Biomed. Opt. Express 1(1), 143–156 (2010). [CrossRef]
24. F. Nouizi, R. Chabrier, M. Torregrossa, and P. Poulet, “3D modeling for solving forward model of no-contact fluorescence diffuse optical tomography method,” Proc. SPIE 7369, 73690C (2009). [CrossRef]
25. S. R. Arridge, “Photon-measurement density functions. Part I: analytical forms,” Appl. Opt. 34(31), 7395–7409 (1995). [CrossRef] [PubMed]
26. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15(2), R41–R93 (1999). [CrossRef]
27. F. Gao, H. Zhao, L. Zhang, Y. Tanikawa, A. Marjono, and Y. Yamada, “A self-normalized, full time-resolved method for fluorescence diffuse optical tomography,” Opt. Express 16(17), 13104–13121 (2008). [CrossRef] [PubMed]
28. W. Becker, The bh TCSPC Handbook, 3rd ed. (Becker & Hickl GmbH, 2008).
29. G. T. Herman, Image Reconstruction from Projections: The Fundamentals of Computerized Tomography (Academic, 1980).
30. A. H. Gandjbakhche, R. F. Bonner, R. Nossal, and G. H. Weiss, “Absorptivity contrast in transillumination imaging of tissue abnormalities,” Appl. Opt. 35(10), 1767–1774 (1996). [CrossRef] [PubMed]
31. A. H. Gandjbakhche, V. Chernomordik, J. C. Hebden, and R. Nossal, “Time-dependent contrast functions for quantitative imaging in time-resolved transillumination experiments,” Appl. Opt. 37(10), 1973–1981 (1998). [CrossRef]
32. R. Cubeddu, A. Pifferi, P. Taroni, A. Torricelli, and G. Valentini, “Imaging of optical inhomogeneities in highly diffusive media: Discrimination between scattering and absorption contributions,” Appl. Phys. Lett. 69(27), 4162–4164 (1996). [CrossRef]
33. R. Cubeddu, A. Pifferi, P. Taroni, A. Torricelli, and G. Valentini, “Imaging with diffusing light: an experimental study of the effect of background optical properties,” Appl. Opt. 37(16), 3564–3573 (1998). [CrossRef]
34. M. S. Patterson, B. Chance, and B. C. Wilson, “Time resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties,” Appl. Opt. 28(12), 2331–2336 (1989). [CrossRef] [PubMed]
35. V. Venugopal, J. Chen, F. Lesage, and X. Intes, “Full-field time-resolved fluorescence tomography of small animals,” Opt. Lett. 35(19), 3189–3191 (2010). [CrossRef] [PubMed]
