Spatially varying regularization based on spectrally resolved fluorescence emission in fluorescence molecular tomography

Johan Axelsson; Jenny Svensson; Stefan Andersson-Engels

doi:10.1364/OE.15.013574

1. Introduction

Fluorescence molecular tomography (FMT) has evolved during the last decade to become an important modality for studying fluorescent markers inside small animal models. Disease and treatment progression can be monitored by imaging the fluorescence emitted from fluorophores attached to specific molecules [1]. One major challenge in FMT is to obtain a feasible tomo-graphic reconstruction resembling the correct result as well as possible. The reconstruction is mathematically ill-conditioned yielding multiple non-unique solutions satisfying the recorded signals. Hence there is a need for methods to increase the robustness of the mathematical problem. This can be achieved by employing a priori information [2, 3]. Prior information is typically some information that imposes a constraint to the problem at hand. Spatial priors can be incorporated in several ways. Anatomical images retrieved from MR has been used to select a feasible reconstruction region in diffuse optical tomography [4]. Further the use of varying regularizing parameters in the reconstruction has been reported to yield better image quality [5].

Recent development of non-contact detection schemes for FMT [6, 7] has provided possibilities for high resolution reconstructions. Using such a scheme the computational cost is greatly increased hence it is important to optimally define the reconstruction parameters balancing computational expense and resolution in the reconstruction[7].

In this paper we report on a simple approach to render a spatially varying regularization map that serves as prior information to the inverse problem of finding the fluorophore absorption coefficient. The method is based on the multispectral fluorescent emission information, as a fluorophore is typically characterized by a broad band fluorescence spectrum. Initial reconstructions of inclusions in a tissue phantom, with and without the prior information adapted, are presented.

2. Theory

2.1. Forward light propagation model

In an FMT setup laser light is radiating several small spots on the tissue surface sequentially, either through a fiber or by a free laser beam. The excitation light propagates through the tissue and undergoes multiple scattering and absorption events. Upon absorption by a fluorophore, fluorescent light will be emitted. The emitted fluorescence then propagates through the tissue and is collected by a detector at the rear boundary. By using multiple detectors one will obtain information on the spatial intensity distribution of the light escaping from the tissue surface.

The steady-state diffusion approximation to the radiative transport equation is utilized to mathematically describe the continuous-wave light propagation in tissue. The light delivered as a small spot onto the tissue surface is then approximated by a point source positioned one scattering length 1/μ′_s inside the medium. The extrapolated boundary condition is applied using mirrored sources to account for the refractive index mismatch at the boundaries [8]. The fluence rate of the excitation light (u_x [W/m ²]) at an arbitrary position r due to a point source at position r _s is then analytically given by the Green’s solution to the homogeneous diffusion equation, Eq. (1).

u_{x} (r_{s}, r) = \frac{P_{0}}{D_{x}} G_{x} (r_{s}, r)

where

G_{x} (r_{s}, r) = \frac{\exp (- μ_{effx} ∣ r - r_{s} ∣)}{4 π ∣ r - r_{s} ∣} .

Here the subscript x denotes excitation wavelength. P ₀ [W] is the laser source power, D_x [m] and μ_effx [m ^-1] is the diffusion coefficient and effective attenuation coefficient, respectively, defined in Eq. (3) and Eq. (4).

D_{x, m} = \frac{1}{3 (μ_{s_{x, m}}^{'} + α μ_{a_{x, m}})}

μ_{{eff}_{x, m}} = \sqrt{\frac{μ_{a_{x, m}}}{D_{x, m}}}

In Eq. (3) and Eq. (4) μ_a[m ^-1] is the absorption coefficient and μ′_s [m ^-1] is the reduced scattering coefficient. As before the subscript x denotes excitation while subscript m denotes emission wavelengths. The dimensionless constant α (here α ≈ 0.55) is adopted from Ripoll et. al. [9] to form a modified absorption dependent diffusion coefficient.

Provided that a fluorophore is present at the position r the excitation light will be absorbed and the fluorophore will emit fluorescence. The strength of the emitted fluorescence is dependent on the absorption coefficient and the fluorescent yield of the fluorophore. The fluorescent yield is characterized by the broad fluorescence spectrum and provides information about the probability that an excitation photon will induce a fluorescent photon at a specific spectral band.

The mathematical formulation of fluorescent light detection relies on a coupled diffusion equation between the excitation light and the emitted fluorescent light [10]. The fluorescence fluence rate (u_m [W/m ²]) at a detector position r _d due to an excitation event in a small volume dV at a position r is then analytically described by Eq. (5) [11].

u_{m} (r_{s}, r, r_{d}) = dV u_{x} (r_{s}, r) \times \frac{μ_{af} (r) γ_{m}}{D_{m}} G_{m} (r, r_{d}) .

The notation follows the convention stated in Eq. (1) where subscript m denotes the emission wavelength. μ_af(r) [m ^-1] is the fluorophore absorption coefficient at the excitation wavelength at an arbitrary position r. γ_m [–] is the fluorescence yield at emission wavelength λ_m. Considering an arbitrary spatial distribution of fluorophore the fluorescence fluence rate at a detector position r _d, due to an excitation source in r _s, is a contribution of all fluorescent volume fractions dV throughout the volume. Hence the fluorescence fluence rate is given by an integral, Eq. (6) [11].

U_{m} (r_{s}, r_{d}) = \int_{V} u_{m} (r_{s}, r, r_{d}) dV .

2.2. Reconstruction of fluorophore inclusion

To reconstruct the spatially varying fluorophore absorption coefficient, i.e. μ_af, the normalized Born approach is applied [12]. Utilizing this approach the fluorescence fluence rate normalized with the excitation wavelength fluence rate at a detector position r _d with the excitation source placed at r _s is given by [12]

U_{nb} (r_{s}, r_{d}) = \frac{U_{m} (r_{s}, r_{d})}{u_{x} (r_{s}, r_{d})} = \frac{1}{u_{x} (r_{s}, r_{d})} \times \int_{V} \frac{u_{x} (r_{s}, r) G_{m} (r, r_{d})}{D_{m}} μ_{af} (r) γ_{m} dV .

The same notation is used in Eq. (7) as above. To solve for the unknown, i.e. μ_af(r) the volume integral in Eq. (7) is discretized into N _voxels where voxel j has a volume of ∆V_j [m ³].

U_{{nb}_{i}} (r_{s_{i}}, r_{d_{i}}) \approx \sum_{j = 1}^{N_{voxels}} w_{i, j} μ_{{af}_{j}} γ_{m}

where

w_{i, j} = {Δ V}_{j} \frac{u_{x} (r_{s_{i}}, r_{j})}{u_{x} (r_{s_{i}}, r_{d_{i}})} \frac{G_{m} (r_{j}, r_{d_{i}})}{D_{m}} .

Here the subscript i denotes the source-detector pair. Equation (8) is a set of linear equations that can be expressed in matrix form, i.e. U _nb = WX. The elements of the weight matrix W is stated in Eq. (9), while

x_{j} = μ_{{af}_{j}} γ_{m}

is the unknown quantity to solve for placed in the column vector X. U _nb is a column vector with the number of elements equal to the number of source-dector pairs [12].

To retrieve the solution to the inverse problem, i.e. X, the weight matrix W should be inverted. Using Tikhonov regularization [13] the regularized solution is given by

X = {(W^{T} W + β I)}^{- 1} W^{T} U_{nb} .

β is a regularization parameter that reduces the ill-condition of the inverse problem. The reconstructed image will depend on the regularization parameter since choosing a too high value of β will produce an image with low contrast and low resolution. On the other hand a too small value of β introduces noise [5, 14]. An empirical way to find the near optimal regularization parameter is to start at a high value and iteratively decrease the value until the error norm decreases to a pre-defined limit [2].

2.3. Spatial priors based on fluorescence emission

The regularization parameter β can be spatially varying since the matrix β I is of the size N _voxels × N _voxels. Hence each voxel can have a voxel-specific β-value. The purpose of using a spatially varying regularization parameter is to increase the spatial resolution in regions where the inclusion most likely is positioned. At the same time regions where the inclusion is not positioned should be smoothed to reduce artifacts. To retrieve the spatially varying regulariza-tion parameter, or regularization map, the fluorescence emission from the inclusion itself can be used.

It has previously been reported that the intensity ratio of the fluorescence emission at two wavelengths, here denoted λ_m1 and λ_m2, from a fluorescent inclusion detected at the boundary of the tissue is dependent on the depth of the inclusion [15, 16]. This effect is due to the difference in bulk tissue attenuation, mainly absorption i.e. μ_{a_m1} = μ_{a_m2}. Following the formalism above, the ratio of the two fluorescence signals induced at a location r is given by

u_{mR} (r, r_{d}) = \frac{u_{m 1} (r_{s}, r, r_{d})}{u_{m 2} (r_{s}, r, r_{d})} = \frac{D_{m 2} γ_{m 1}}{D_{m 1} γ_{m 2}} \frac{G_{m 1} (r, r_{d})}{G_{m 2} (r, r_{d})} .

Equation (12) can be calculated to yield the ratio of the escape probabilities for the two wavelengths emitted from any voxel to each detector. This is performed by utilizing the reciprocity theorem [17, 18], i.e. G_m(r,r _d) = G_m(r _d,r). Calculation of Eq. (12) from detector d yields a map of the ratio throughout the geometry. For a given set of optical properties the value of the ratio, in Eq. (12) is only dependent on the distance from the detector. Due to the spherical symmetry of the Green’s solution all voxels positioned at the same distance from the detector will have the same value of the calculated ratio. Comparing the forward modelled ratio in each voxel and the experimentally acquired ratio, i.e. Ũ_mR(r _d), will then result in a minima in all voxels where the concurrence between the modelled ratio and the measured ratio is optimal. The comparison is performed by simply taking the difference between the two quantities which can be expressed as

Δ U_{mR}^{d} (r_{j}) = ∣ In (u_{mR} (r_{d}, r_{j})) - In ({\tilde{U}}_{mR} (r_{d})) ∣ .

Here subscript j denotes voxel number in the discretized geometry and Ũ_mR(r _d) represents the ratio of the measured fluorescence intensity at the two wavelengths. In the evaluation of the difference between modelled and measured values, we have in Eq. (13) chosen to use logarithmic values in order to decrease the dynamic of the resulting difference. Evaluation of Eq. (13) for one detector and all voxels in a plane through the volume yields an arc with minimum values inside the geometry, shown in Fig. 3(a-c). The minima distribution shaped like an arc results from the spherical symmetry in the calculation of Eq. (12). Due to the actual propagation distance from the fluorophore to each detector the measured intensity ratio will be different. A detector collecting fluorescence that have propagated a long distance will yield a small intensity ratio and the arc will be placed far from the detector. Analogously a detector close to the fluorophore will detect a large ratio hence resulting in an arc with smaller radius.

To retrieve a regularization map to be used as a spatially varying regularization parameter the mean of Eq. (13) for all detectors is calculated through

\bar{Δ U_{mR}} (r_{j}) = \frac{1}{N_{detectors}} \sum_{d = 1}^{N_{\det ectors}} Δ U_{mR}^{d} (r_{j}) .

Equation (14) will yield a minimum where the arcs overlap, i.e. where the origin of fluorescent emission is positioned. The intersection of the arcs is shown schematically in Fig. 3(d). Equation (14) is then subject to scaling. The empirically determined scaling method utilized here is stated in Eq. (15).

p_{j, j} = In (\frac{\bar{Δ U_{mR}} (r_{j})}{min (\bar{Δ U_{mR}} (r_{j}))}) + 1 .

where p_j,j is the matrix element of voxel j in the diagonal matrix P. The incorporation of the regularization map P into the reconstruction algorithm is governed by replacing the matrix β I in Eq. (11) with the matrix β P. The regularization parameter will then have a minimum value of β which is iteratively decreased during the reconstruction. The initial value of the parameter is ten times the maximum diagonal element of the matrix W ^T W.β is decreased until the projection error, i.e. |Ũ_nb - WX| /N _voxels, is lower than 0.05.

The number of sources to be used is arbitrary. Here we have used the fluorescence measurements due to a source positioned in front of the inclusion, which yields the highest level of fluorescence signal.

3. Experimental setup

3.1. Imaging system

The experimental setup can be seen in Figure 1. The setup consisits of a fiber-coupled light source, tissue phantom and a multispectral imaging system. The light source was a laser at 532 nm (Millennia Vs, Spectra Physics Laser) emitting at 1 W. The laser light was collimated and coupled into an optical fiber with a core diameter of 400 μm, the power at the distal end of the fiber was approximately 80 mW. This end of the fiber was translated over one of the boundaries in a source grid made of black delrin plastic with drilled holes providing a source position separation of 5 mm. Emitted fluorescence from the object was detected by a multispectral imaging system consisting of a CCD-camera (C4742-80-12AG, Hamamatsu), a liquid crystal tunable filter (LCTF VIS 20–35, Varispec) and a camera objective lens (Nikon f/1.8, focal length 50 mm). The LCTF filter was scanned in the wavelength range of 532–660 nm in steps of 10 nm to collect light. Each spectral band has a spectral width of 20 nm.

Fig. 1. The experimental setup
Download Full Size | PDF

3.2. Tissue phantom

The tissue-simulating phantom consisted of a glass cuvette with a thickness of 2 cm, well in the range of small animal imaging studies. The solution with tissue like optical properties was a solution of Intralipid (Fresenius Kabi, Sweden), water and bovine blood. The scattering is due to the Intralipid with a concentration of 1.1%. The absorption of the phantom is due to the bovine blood. The absorption coefficient and reduced scattering coefficient of the two phantom solutions were evaluated using an integrating sphere setup [19]. The absorption coefficient and reduced scattering coefficient can be seen in Fig. 2(a) and Fig. 2(b) respectively. A glass cylinder with inner diameter (d_cyl=5.6 mm) was filled with the phantom solution and Rhodamine 6G, a fluorescent dye with similar emission spectra to fluorescent proteins. Two concentrations of Rhodamine 6G were used c ₁=0.5 μM and c ₂=1 μM. The emission characteristics of the fluorophore was acquired using a spectrometer (USB2000, OceanOptics) and the normalized spectrum is shown in Fig. 2(c).

Fig. 2. Optical properties for the experimental tissue phantom obtained from integrating sphere measurements. In (a) absorption coefficient and (b) the reduced scattering coefficient. In (c) the normalized fluorescence spectrum for the fluorophore is shown.
Download Full Size | PDF

3.3. Measurement procedure and data analysis

The glass cylinder was placed inside the cuvette at a depth of z = 7 or z = 11 mm. Images were acquired with the filter centered at the excitation wavelength 532 nm and at the following fluorescence wavelengths 550–660 nm in steps of 10 nm. Background images were acquired with no excitation light. Images with no fluorescent object inserted into the phantom were also obtained to characterize the autofluorescence from the phantom material. For all measurements the exposure times were in the range of 2–6 seconds.

All images were subtracted with the corresponding background image and normalized with the exposure time. The background fluorescence from the phantom material was also subtracted in order to provide data free from autofluorescence. A total number of 10 source positions were used and 31 virtual detectors were extracted from the images after binning of 4×4 pixels. In the data processing level detected data was normalized with the spectrally varying CCD quantum efficieny and filter transmission. Measurements with lower signal strength than the background level was rejected.

The fluorescent yield was retrieved prior to the data evaluation using the emission spectrum in Fig. 2(c). Due to the spectral width of the filters used the fluorescent yield was calculated using

γ_{m} = Q_{e} \frac{\int_{λ_{m} - 10}^{λ_{m} + 10} γ (λ) dλ}{\int_{0}^{\infty} γ (λ) dλ}

where Q_e = 0.95 is the quatum yield for Rhodamine 6G [20] and γ(λ) is the normalized fluorescence spectrum seen in Fig. 2(c).

4. Results and Discussion

4.1. Spatially varying regularization map based on fluorescence emission

The cylinder with a concentration of 1 μM was placed at z = 11 mm. Calculation of Eq. (13) using emission measurements at λ_m1 = 560 nm and λ_m2 = 600 nm for three detectors yields the cross-sectional images shown in Fig. 3. The measurements were extracted with the source position fixed, indicated in Fig. 3. The resulting arc is centered around the specific detector, since the minimum is obtained for all voxels positioned at the same distance from the detector. The minimum value for the three detectors shown in the figure are placed at different distances from the detectors due to the fact that the measured ratios Ũ_mR are different, requiring a matching value of |r _j -r _d| in the forward model u_mR to compensate. A longer propagation distance yields a smaller intensity ratio, hence the arc is placed further away from the detector.

Fig. 3. ∆U_mR for detector (a) 5, (b) 15 and (c) 27 extracted for the same excitation source position. The source position (•) and detector position (x) is indicated. The true position of the cylinder is marked by the circular ring for reference. In (d) the summation of the minima arcs in (a)-(c) are shown schematically. The colouring indicate the detector each arc is originating from.
Download Full Size | PDF

The spatially varying regularization maps retrieved using Eq. (14) for four different cases is shown in Fig. 4. The maps rendered from the cylinder with 0.5 μM concentration placed at z = 11 mm z = 7 mm is shown in Fig. 4(a) and (b), respectively. The regularization maps for the case with 0.5 μM concentration placed at z = 11 mm z = 7 mm is shown in Fig. 4(c) and (d), respectively. All regularization maps retreived through the emission ratio within this paper are based on λ _m1 = 560 nm and λ_m1 = 600 nm. The higher absorption for λ _m1 will produce maps with higher spatial resolution as seen in Fig. 4(b) and (d). Inclusions positioned far from the detection boundary will also have a slightly more narrow intensity distribution over the boundary. This is seen in the maps where Fig. 4(b) and (d) are extended in x-direction while Fig. 4(a) and (c) are more symmetric around the inclusion. The method utilizes a ratio of two emission spectral bands where the dependence on the fluorophore concentration disappears. This is also seen in Fig. 4(a)–(d) where the regularization parameter values are effectively the same. Small differences are due to experimental fluctuations.

Fig. 4. Regularization maps retrieved using the ratio of two spectral bands within the fluorescence emission profile. The fluorophore concentration (c) and position (z) is (a) c = 0.5 μM and z= 11 mm, (b) c =0.5 μM and z = 7 mm, (c) c = 1 μM and z=11 mm, (d) c=1 μM and z = 7 mm. A circular ring indicates the true inclusion position.
Download Full Size | PDF

4.2. Reconstruction of single inclusion

The reconstruction with no prior information of a cylinder with concentration 1 μM positioned at z = 11 mm is shown in Fig. 5(a). The same reconstruction using the spatially varying regu-larization map in Fig 4(c) is shown in Fig. 5(b). The cross-sectional plots in z- and x-direction is inserted in Fig. 5(c) and (d). The target fluorophore absorption coefficient was μ_af = 23 m ^-1 indicated in Fig. 5(c) and (d). It is seen that the prior regularization map reduces the artifacts around the inclusion as well as around source and detector boundaries. This effect is due to the higher value of the regularization parameter further away from the inclusion which effectively smoothes the solution.

The cross-sectional plots in Fig. 5(c) and (d) shows that the reconstructed fluorophore absorption coefficient is closer to the true value using prior regularization map. It is also seen that the cross-section is more narrow while the reconstruction without prior yields a broad inclusion cross-section. The method presented herein scales the regularization map and relates all values in the map to one voxel, i.e. the minimum of the maps shown in Fig. 4. The regularized solution will then have a maximum in that voxel and the spatial resolution of the reconstruction will depend on the scaling of the regularization parameters. Fixing the regularization parameter on the boundaries and decreasing the minimum regularization parameter value, through scaling, Pogue et. al. [5] reported that a higher contrast will be retrieved in diffuse optical tomographic reconstructions. This can be achieved using a modification of the scaling in Eq. 15.

Fig. 5. Reconstructions of a cylinder with a concentration of 1 μM positioned at z = 11 mm using (a) no prior and (b) spatially varying regularization map. In (c) and (d) the cross-sections indicated by dotted lines in (a) and (b) are shown. The full line (–) represents the reconstruction with prior information while the dashed line (––) represent reconstruction without prior. A circular ring indicates the true position of the inclusion.
Download Full Size | PDF

4.3. Reconstruction of two inclusions

Two cylinders of concentrations 0.5 and 1 μM were placed at a depth of z = 9 mm and imaged using the imaging system. To render the regularization map several source positions were used. Equation (15) was then calculated twice; one with the fluorescence measurements acquired using source positions in the interval x_source ∈ [0,25] yielding P ₁ and the other with x_source ∈ [30,50] resulting in P ₂. The selection of source positions are chosen based on the intensity distribution of the fluorescence emission at the detector side. Hence the regularization map P ₁ is based on only those measurement where one of the fluorophores is excited while the excitation of the other inclsuion result in P ₂. The two regularization maps are then added together by selecting the minimum value of either P ₁ or P ₂ for a specific voxel. The resulting regularization map is seen in Fig. 6(b). The reconstructions were performed as above and the result without prior information is seen in Fig. 6(a) and with prior information in Fig. 6(c). As before the artifacts are reduced when reconstructing with prior information. Quantitatively the prior has less impact. This is due to the small dynamics of the scaling ranging only from ~ 1.5 to ~ 3.5, see Fig. 6(b). The reconstructed values of the two inclusions are in both cases, i.e. with and without prior, about 50 % of the true value.

5. Conclusion

We have presented a simple and objective procedure that renders prior information about the position of a fluorescent inclusion based on the multispectral feature of the fluorescence. The method is based on the difference of the forward model and the detected intensities. This poses no additional computational cost for a reconstruction algorithm, except the forward matrix calculation of the additional emission wavelength. Based on initial image reconstructions the conclusion is that the most prominent feature of using the spatially varying regularization parameter is the reduction of the artifacts close to source and detector boundaries.

Fig. 6. Reconstructions of two cylinder with concentrations 1 μM (lower) and 0.5 μM (upper) positioned at z = 9 mm using (a) no prior and (c) spatially varying regularization map. In b) the regularization map is shown. The circular rings indicate the true position of the inclusions.
Download Full Size | PDF

The scaling of the regularization map affects the solution. The reported method finds the minimum in the regularization map, outlining the most probable inclusion position, and scales all other voxel values with respect to the minimum. The scaling equation used here is empirically retrieved and in future work we will investigate the use of modifications to this equation in order to find an optimal scaling method.

The demonstration presented herein is performed in a homogeneous tissue phantom with autofluorescence subtracted from the measurements. The influence of autofluorescence and heterogeneous media will most likely affect the retrieval of the regularization maps since these are based on the emitted fluorescence detected at the boundary. Heterogeneous media increases the complexity of the inverse problem hence it is of even greater importance to render prior information under these circumstances. In what respect these factors influence the regularization map is now being investigated in our present work.

Acknowledgements

We acknowledge the discussions with Oliver Dorn. This research was supported by the EU Integrated Project Molecular Imaging LSHG-CT-2003-503259.

References and links

1. V. Ntziachristos, J. Ripoll, L. Wang, and R. Weissleder, “Looking and listening to light: the evolution of whole-body photonic imaging,” Nat. Biotechnol. 23, 313–320 (2005). [CrossRef] [PubMed]

2. B. Brooksby, S. Jiang, H. Dehghani, B. Pogue, K. Paulsen, J. Weaver, C. Kogel, and S. Poplack, “Combining near-infrared tomography resonance imaging to study in vivo and magnetic breast tissue: implementation of a Laplacian-type regularization to incorporate magnetic resonance structure,” J. Biomed. Opt. 10 (2005). [CrossRef] [PubMed]

3. M. Guven, B. Yazici, X. Intes, and B. Chance, “Diffuse optical tomography with a priori anatomical information,” Phys. Med. Biol. 50, 2837–2858 (2005). [CrossRef] [PubMed]

4. H. Xu, R. Springett, H. Dehghani, B. Pogue, K. Paulsen, and J. Dunn, “Magnetic-resonance-imaging-coupled broadband near-infrared tomography system for small animal brain studies,” Appl. Opt. 44, 2177–2188 (2005). [CrossRef] [PubMed]

5. B. Pogue, T. McBride, J. Prewitt, U. Osterberg, and K. Paulsen, “Spatially variant regularization improves diffuse optical tomography,” Appl. Opt. 38, 2950–2961 (1999). [CrossRef]

6. R. Schulz, J. Ripoll, and V. Ntziachristos, “Experimental Fluorescence Tomography of Tissues With Noncontact Measurements,” IEEE Trans. Med. Imaging 23, 492–500 (2004). [CrossRef] [PubMed]

7. E. Graves, J. Culver, J. Ripoll, R. Weissleder, and V. Ntziachristos, “Singular-value analysis and optimization of experimental parameters in fluorescence molecular tomograpy,” J. Opt. Soc. Am. 21, 231–241 (2004). [CrossRef]

8. D. Contini, F. Martelli, and G. Zaccanti, “Photon migration through a turbid slab described by a model based on diffusion approximation . 1. Theory,” Appl. Opt. 36, 4587–4599 (1997). [CrossRef] [PubMed]

9. J. Ripoll, D. Yessayan, G. Zacharakis, and V. Ntziachristos, “Experimental determination of photon propagation in highly absorbing and scattering media,” J. Opt. Soc. Am. 22, 546–551 (2005). [CrossRef]

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

11. M. O’Leary, D. Boas, X. Li, B. Chance, and A. Yodh, “Fluorescence lifetime imaging in turbid media,” Opt. Lett. 21, 158–160 (1996). [CrossRef] [PubMed]

12. V. Ntziachristos and R. Weissleder, “Experimental three-dimensional fluorescence reconstruction of diffuse media by use of a normalized Born approximation,” Opt. Lett. 26, 893–895 (2001). [CrossRef]

13. A. Neumaier, “Solving ill-conditioned and singular linear systems: A tutorial on regularization,” Siam Review 40, 636–666 (1998). [CrossRef]

14. H. Dehghani, D. Barber, and I. Basarab-Horwath, “Incorporating a priori anatomical information into image reconstruction in electrical impedance tomography,” Physiol. Measurement 20, 87–102 (1999). [CrossRef]

15. J. Swartling, J. Svensson, D. Bengtsson, K. Terike, and S. Andersson-Engels, “Fluorescence spectra provide information on the depth of fluorescent lesions in tissue,” Appl. Opt. 44, 1934–1941 (2005). [CrossRef] [PubMed]

16. J. Svensson and S. Andersson-Engels, “Modeling of spectral changes for depth localization of fluorescent inclusion,” Opt. Express 13, 4263–4274 (2005). [CrossRef] [PubMed]

17. S. Arridge, “Optical tomography in medical imaging,” Inverse Problems 15, R41–R93 (1999). [CrossRef]

18. J. Swartling, A. Pifferi, A. Enejder, and S. Andersson-Engels, “Accelerated Monte Carlo model to simulate fluorescence spectra from layered tissues,” J. Opt. Soc.Am. 20, 714–727 (2003). [CrossRef]

19. J. Dam, T. Dalgaard, P. Fabricius, and S. Andersson-Engels, “Multiple polynomial regression method for determination of biomedical optical properties from integrating sphere measurements,” Appl. Opt. 39, 1202–1209 (2000). [CrossRef]

20. R. Kubin and A. Fletcher, “Fluorescence Quantum Yields of Some Rhodamine Dyes,” J. Lumin. 27, 455–462 (1982). [CrossRef]