## Abstract

Experimental results for imaging the low-scattering tissue phantoms based on the derivative estimation through perturbation Monte-Carlo (pMC) method are presented. It is proven that pMC-based methods give superior reconstructions compared to diffusion-based reconstruction methods. An easy way to estimate the Jacobian using analytical expression obtained from perturbation Monte-Carlo method is employed. Simulation studies on the same objects, considered in the experiment, are performed and corresponding results are found to be in reasonable agreement with the experimental studies. It is shown that inter-parameter cross talk in diffusion based methods lead to false results for the low-scattering tissue, where as the pMC-based method gives accurate results.

©2005 Optical Society of America

## 1. Introduction

Biomedical optical imaging [1–8] in the recent years has been the topic of much interest to many researchers around the globe. Near Infra-red (NIR) optical imaging is mainly used to image thick tissue such as human brain, breast and joints. Main aim of NIR imaging is to get a spatial distribution of optical properties from the measurements done on the boundary of tissue. Light propagation in tissue can be accurately described by Radiative Transfer Equation (RTE) [9], which is an integro-differential equation. RTE is difficult to solve either analytically or numerically. The diffusion approximation [10] is generally applied to RTE and solved therein to reconstruct the optical properties. Even though diffusion approximation will be accurately able to describe the light propagation in tissue in most of the cases, it fails in the cases of low-scattering tissue [11–14], such as Cerebral Spinal Fluid (CSF) layer of the brain [15–19] and Synovial fluid layer in the joints [20–22]. Numerical solutions from RTE are proven to be computationally expensive [23–27]. Hybrid radiosity finite-element- and higher-order diffusion models [28–30] require *a priori* knowledge of the boundaries of the low-scattering regions. This *a priori* knowledge is difficult to obtain in most of the practical situations, for example neonatal brain grows immensely in both shape and complexity immediately after birth [31]. One of the main issues of diffuse optical tomography is the inter parameter cross-talk [32–34], which is due to underlying non-uniqueness in the inverse problem.

To image these low scattering tissues, there have been attempts to incorporate Monte Carlo (MC) simulation for doing reconstruction problem [35–37] (*inverse problem* of optical tomography). MC simulation is considered to be equivalent of implementing the RTE, which is accurate to describe the light propagation in these cases. Here again, MC simulation takes large computational time involved for taking many millions of photons through the tissue. Inversion methods that make use of repeated application of the forward model and updating of the derivative needed for guiding the solution to convergence would become too expensive to be a useful option. This problem has been recently addressed [36,37] making use of the perturbation Monte-Carlo (pMC) method to extract Jacobian, which is the derivative of measurements (forward model) with respect to the optical absorption (µ_{a}) and scattering coefficient (µ_{s}).

An analytical expression for the perturbation in the detected photon weight consequent to changes in optical properties in a sub-region, borrowed from neutron-transport theory, is used to update both Jacobian and the computed forward data [36,37]. After discretization, such a procedure will require only a single Monte-Carlo (MC) simulation with the derivative- as well as forward projection data extraction handled by the pMC, which is rapidly done. In Ref. [36], this is used to solve a simple two region inverse problem of photon migration in a heterogeneous slab of thickness comparable to the transport mean free path, an object which falls under the transport-regime, where diffusion equation (DE) fails to hold. An extension of the pMC approach to construct a Jacobian matrix with basic MC simulation as forward model is presented in Ref. [37] for use in a perturbation-based method to reconstruct transport-regime low-scattering objects with discretized absorption- and scattering coefficient distributions.

With numerical simulations the utility of a pMC-based iterative reconstruction strategy is shown [37] to reconstruct low-scattering objects with single- or multiple inhomogeneities. Verification for the utility of the above algorithm by reconstructing µ_{a}- and/or µ_{s} distributions from experimental time-domain data obtained from low-scattering phantoms with single- or multiple inhomogeneities is presented here. The propagation-backpropagation (PBP) strategy [38], where data is handled one source location at a time is used in the iterative reconstruction scheme, which reduces the dimension of the problem and the overall computation time. Reconstructed images of tissue-equivalent phantoms using the experimental data are presented. Even with the total-intensity measurements it is demonstrated that, one can clearly distinguish absorption- and scattering heterogeneities, in which case, diffusion based reconstruction methods require additional constrains [33].

The work presented here is as follows:- In section 2, a brief overview of the iterative reconstruction scheme is presented which uses the Jacobian to calculate the update vectors. A brief description of the calculation of the Jacobian matrix using the pMC is also given. Then, in section 3, the details of the time-domain experiments conducted are described. The tissue-equivalent phantoms and the geometry of data collection are also given in section 3. In section 4, reconstruction results from the experimental data along with simulation results are presented. Finally, to show the effectiveness of this approach, one set of reconstruction results from the simulated data using the diffusion-equation (DE) based reconstruction method [39] are also presented.

Here only PBP method is used where data from a single source (and many detectors) are reconstructed to arrive at an updated object which is used as the initial estimate for reconstructing data from the next source. In Ref. [37], the effectiveness of the PBP approach is shown for doing inverse problem compared to conventional approach. Since the number of equations in the PBP approach is smaller, the reconstruction problem will be too ill-posed to converge to a solution when both scattering- and absorption inhomogeneities are considered simultaneously without *a priori* knowledge on location [37]. For reduction in ill-posedness of the inverse problem and for better convergence, an *a priori* information is used in all the cases, except one, presented here. Simulation studies for these objects are also presented in this work (details are given in Section 4). As a part of the simulation studies, it is shown that with or without *a priori* information, a pMC-based approach will be able to give same quality of reconstructions.

## 2. Iterative reconstruction algorithm

Complete details of the iterative reconstruction algorithm are presented in Ref. [37]. Here this is briefly explained. The steps involved in the iterative reconstruction procedure are shown in the flow chart of Fig. 1. It has two iterations, one main one (the outer iteration) and another subsidiary one (the inner iteration) [40]. The i^{th} solution vector (${\mu}_{\mathrm{a}}^{\mathrm{i}}$, ${\mu}_{\mathrm{s}}^{\mathrm{i}}$) given to the forward model (the Monte-Carlo procedure is repeated with 2.1 million photons to model the forward propagation) predicts the computed data (W^{c}). The experimental data (W^{e}) plugged in help us find ΔW(=W^{e}- W^{c}) that is used to set-up the update-equation (Eq. (1)). In this work, integrated (or DC) intensity (Eq. (4.8) in Ref. [8]) is used as the data. The inner iteration outputs the update vector (Δµ_{a}, Δµ_{s}), which is used to arrive at the new-solution ${\mu}_{\mathrm{a}}^{\mathrm{i}+1}$=${\mu}_{\mathrm{a}}^{\mathrm{i}}$+Δµ_{a} and ${\mu}_{\mathrm{s}}^{\mathrm{i}+1}$=${\mu}_{\mathrm{s}}^{\mathrm{i}}$+Δµ_{s}. The conjugate gradients squared (CGS) method is used to solve this minimization problem (inner iteration).

Analytic expressions for changes introduced in detected photon weights owing to optical property perturbation in a sub-region in the object are made use of to calculate both the rate of change of photon weights with respect to optical properties and the new photon weights [36,37,41]. Specifically if µ_{a} and µ_{s} are perturbed in certain locations and become *$\overline{\mu}$ _{a}*=

*µ*

_{a}+

*δµ*

_{a}and

*$\overline{\mu}$ =µ*then the detected photon weight

_{s}+δµ_{s}*w*changes to

*w̄*as

where *n* is the number of collisions the photon undergoes in the perturbed region, and *l* is the path-length traversed therein and µ_{t}=µ_{a}+µ_{s}. Eq. (1) provides a way to estimate the changes in measured photon density (or weight) owing to changes in µ_{a} and µ_{s} in a specified location.

Here the object is discretized into 81×81 pixels with freedom for µ_{a} and µ_{s} values in these pixels to move towards their values at convergence. If data from only M detectors is considered for each of the S source positions the Jacobian matrix connecting the change in measurements at the boundary to perturbations in µ_{a} and µ_{s}, has dimensions either (S×M)·{(81×81)×2} or M×{(81×81)×2} (when data from each source is handled separately which is the PBP approach), depending on whether the handling of the data is from all the sources together or one source at a time. For constructing such a Jacobian, using Eq. (1), there are 2×(81×81) regions of possibly different µ_{a} and µ_{s} values centered around each pixel. Approximately circular sub-regions containing 12 pixels are introduced surrounding each pixel in the original object domain (for the boundary pixels this is done by extending the boundary), with the object assigned homogeneous background values, ${\mu}_{\mathrm{a}}^{\mathrm{b}}$ and ${\mu}_{\mathrm{s}}^{\mathrm{b}}$ for the calculations of the Jacobian. A large number of photons (28.8 million) are taken through the object to determine the average number of collisions, *n*, and the length traversed, *l*, in each of the sub-regions entered around each of the pixels for the detected photons. The set of *n* & *l* values determined are frozen and used to update the Jacobian matrix during the course of the iterations. To arrive at a particular element of the Jacobian matrix, which is the rate of change of data at a detector with respect to the optical properties at a certain pixel, the stored *n* & *l* values corresponding to that pixel sub-region are used to determine $\frac{\partial \stackrel{\u0305}{w}}{\partial \delta {\mu}_{a}}$ and $\frac{\partial \stackrel{\u0305}{w}}{\partial \delta {\mu}_{s}}$ using Eq. (1). More details about the weight matrix (Jacobian) generation can be found in Section 3 of Ref. [37].

## 3. Experimental details

#### 3.1 Phantom fabrication

Tissue-equivalent phantoms are fabricated as described below: The background material is obtained by mixing 100 parts of Lapox A-53 epoxy resin with 10 parts of K-6 hardener (Cibatul Limited, Gujarat, India). Scattering is introduced by adding TiO_{2} powder (Dupont, Ti-Pure R-902) and absorption through India ink of appropriate quantities as required.

The phantoms with inhomogeneities are fabricated in two steps: In the first a cylinder is cast which is designed to have background optical properties, ${\mu}_{\mathrm{a}}^{\mathrm{b}}$=0.008 mm^{-1} and ${\mu}_{\mathrm{s}}^{\mathrm{b}}$=0.05 mm^{-1}. For this, while mixing the epoxy resin to hardener 469 mg of TiO_{2} powder and 8 mg of 2% India ink are also added and thoroughly mixed. The resin is allowed to set for 24 hours at room temperature (approx. 26°C). Inhomogeneity is introduced as smaller diameter rods of different optical properties. In the experiments conducted here, the inhomogeneity values at its centre in absorption (${\mu}_{\mathrm{a}}^{\text{in}}$) and scattering (${\mu}_{\mathrm{s}}^{\text{in}}$) are varied from 0.008–0.021 mm^{-1} and 0.005–0.14 mm^{-1} respectively. For this in a 110ml mixture of epoxy resin and hardener, 262mg TiO_{2} and 1.6mg 2% Indian ink are added. Both the background cylinder and the rods are thoroughly degassed in a vacuum chamber. The rod, which is to serve as the inhomogeneity, is machined and fitted exactly into a hole drilled in the background cylinder. A composite object with multiple inhomogeneities is also prepared following the same procedure.

The details of tissue-equivalent phantoms fabricated using the above procedure are as follows: Overall diameter=80 mm. Inclusion diameter=10 mm. The background optical properties, ${\mu}_{\mathrm{a}}^{\mathrm{b}}$=0.008 mm^{-1} and ${\mu}_{\mathrm{s}}^{\mathrm{b}}$=0.05 mm^{-1}. The optical properties of the absorption-and scattering inclusion are ${\mu}_{\mathrm{a}}^{\text{in}}$=0.021 mm^{-1} and ${\mu}_{\mathrm{s}}^{\text{in}}$=0.14 mm^{-1} respectively. The refractive index of the material, which is considered uniform, is 1.53. Henyey-Greenstien phase function is used to describe the scattering. The anisotropy factor g is kept at 0.75.

#### 3.2 Data gathering

A schematic of the setup used for experiments is shown in Fig. 2. The second harmonic output of a ps Nd:YAG laser (Continuum PY61C-10, 532 nm, 27 ps, 10 Hz) is used to pump a 100 cm long Raman cell filled with H_{2} at 30 Atm. pressure to generate 683 nm Stokes output by Stimulated Raman Scattering (SRS). Typical conversion efficiency of the set up is ~20% with pulse-to-pulse fluctuation of ~8% in first stokes output at 683nm. The Stokes output is separated from pump and the higher order Stokes and anti-Stokes output by proper high and low pass filters. The spatially filtered Stokes beam is used to illuminate the object mounted on a jig with provision for rotating to any angle. A multimode fiber (*Φ*=1000 µm) of 1 meter length is used for collecting the scattered light exiting from the phantom. The detector fiber position can be independently adjusted to any position around the phantom. The Stokes light emerging from the scattering medium gets temporally elongated due to multiple scattering. A single shot streak camera [42] is used to record the photon arrival histogram (the TPSF). Data is collected for 13 detector positions, with an angular spacing of 5° covering an angle of 60° diametrically opposite to the source position for each of 12 source positions at spacing of 30° covering the phantom fully. Averaging of 50 data sets for each source-detector location is performed, after making suitable correction for temporal jitter, before the data is subjected to further analysis.

## 4. Experimental and simulation results

#### 4.1 Experiments

An initial Monte-Carlo simulation with 28.8 million photons having an extended boundary and homogeneous background values of ${\mu}_{\mathrm{a}}^{\mathrm{b}}$ and ${\mu}_{\mathrm{s}}^{\mathrm{b}}$ is performed to store *n* & *l*, the average number of collisions and length of traverse of photons, in the 12-pixels circular sub-regions surrounding each of the pixels. This data set is needed for calculation of the Jacobian as given in section 2 (see Fig. 1). For forward computed data (W^{c}), integrated intensity, an MC simulation with 2.1 million photons is done on the object with same source-detector positions with the same detector diameter. PBP strategy is used in combination with *a priori* information to reduce the computation effort. Here (in PBP) data from one view is considered sufficient to solve the iteration, which resulted in an updated object for use with data from the next view. All the views are considered cyclically. With µ_{a}- and µ_{s} inhomogeneities simultaneously present, without assumption of locations, the PBP-based problem is too ill-posed to result in a solution [37]. The convergence criterion used is that the norm of the difference between computed data and the experimental data should be below a, preset, small value. Since Monte-Carlo simulation is used, a stochastic forward model, the reconstructed images required post-processing. A local median filter of dimension 5×5 pixels is employed to smoothen the images.

As stated earlier, experimental data (integrated intensity) is collected for three phantoms: first one with only absorption inhomogeneity, the second only with scattering inhomogeneity and third one with both at different locations. In each of the cases, the TPSF’s measured are integrated to obtain the total photon weight (integrated intensity) measured for that source-detector pair (to serve as W^{e}). Data from a particular view is input to the iteration of Fig. 1, which outputs the update vector [Δµ_{a}, Δµ_{s}]. The updated µ_{a} and µ_{s} are used along with the data from the next ‘view’ to continue the iterations (PBP approach). In all the experimental cases considered, it is assumed that the location of inhomogeneity is *a priori* known to be within a region of 30×30 pixels. The initial estimates of (µ_{a}, µ_{s}) for the reconstruction algorithm are the background values (${\mu}_{\mathrm{a}}^{\mathrm{b}}$, ${\mu}_{\mathrm{s}}^{\mathrm{b}}$). All the reconstruction procedures are carried out on a Pentium IV 2.40 GHz processor.

Figure 3 gives the reconstruction of µ_{a}-distribution from the experimental data obtained from the tissue phantom with only µ_{a}-inhomogeneity. The unknowns are the value of µ_{a} at these spatial *a priori* known locations. What is assumed *a priori* is that the inhomogeneity is somewhere inside a region of 30×30 pixels centered at (50, 41) making number of unknowns as 900. As PBP approach is used Jacobian for this problem has a dimension of 13×900. The original object shown in Fig. 3(a) has a µ_{a}-inhomogeneity. The reconstructed µ_{a}-distribution is as shown in Fig. 3(b). The solution converged in 43 iterations, each iteration taking ≈16 minutes. The original object µ_{a}-inhomogeneity value is 0.021 mm^{-1} and the reconstructed ${\mu}_{\mathrm{a}}^{\text{in}}$ at the center is 0.028 mm^{-1}.

Figure 4 shows the result for a similar effort, as above, done on experimental data from the phantom with only µ_{s}-inhomogeneity. Original µ_{s}-inhomogeneity value is 0.14 mm^{-1} (Fig. 4(a)) and reconstructed ${\mu}_{\mathrm{s}}^{\text{in}}$ at the center is 0.28 mm^{-1} (Fig. 4(b)). The solution converged in 59 iterations, each iteration taking ≈16 minutes.

In Fig. 5, both µ_{a}- and µ_{s} inhomogeneities are located at different *a priori* known locations (as before, the locations specified to be within 30×30 pixels regions, for µ_{a} centered at (40, 41) and for µ_{s} at (50, 41)). Here ${\mu}_{\mathrm{a}}^{\mathrm{b}}$ is assumed to be 0.008 mm^{-1} and ${\mu}_{\mathrm{s}}^{\mathrm{b}}$=0.05 mm^{-1}. Note that only PBP approach is used for reconstruction. There are 900×2 unknowns and the Jacobian size is 13×1800. The original ${\mu}_{\mathrm{a}}^{\text{in}}$(=0.021 mm^{-1}) and ${\mu}_{\mathrm{s}}^{\text{in}}$(=0.14 mm^{-1}) images are shown in Fig. 5(a) & (b) respectively. The reconstructions are shown in Fig. 5(c) & (d) with the central value of ${\mu}_{\mathrm{a}}^{\text{in}}$- and ${\mu}_{\mathrm{s}}^{\text{in}}$ reconstructions 0.028 mm^{-1} and 0.21 mm^{-1} respectively. The solution converged in 50 iterations, each iteration taking ≈32 minutes.

#### 4.2 Simulations

Results of the numerical simulations for the same objects considered to those used in the experiments are presented here. The detectors are considered to be of diameter 1mm and simulated experimental data, (W^{e}), are generated by adding 2% Gaussian noise to the integrated intensity arrived at using MC simulation. Similar to the experimental reconstruction procedure, the homogeneous background optical properties (${\mu}_{\mathrm{a}}^{\mathrm{b}}$, ${\mu}_{\mathrm{s}}^{\mathrm{b}}$) are considered as the initial estimate to start the iteration and Pentium IV 2.40 GHz processor for computations.

In the numerical simulation, for the object with only absorption inhomogeneity (Fig. 3(a)), using the PBP approach and the *a priori* assumption of inhomogeneity location, convergence is obtained in 16 iterations. The result is shown in Fig. 6(a). Without *a priori* information, the reconstruction problem took longer time (29 iterations, each taking ≈43 minutes) to converge and the result, shown in Fig. 6(b), is comparable in accuracy to the one shown in Fig. 6(a). Reconstructed µ_{a}-inhomogeneity value at the center of Fig. 6(a) and 6(b) are 0.018 mm^{-1} and 0.019 mm^{-1} respectively. For this reason, in all the other cases considered here, *a priori* location information is assumed to be known. The reconstruction result of simulations corresponding to the object shown in Fig. 4(a), is given in Fig. 7. Convergence is obtained in 21 iterations, with the µ_{s}-value at the center of reconstructed inhomogeneity 0.16 mm^{-1}. Figs. 8(a) and 8(b) show the simulation results for the objects of Fig. 5(a) and 5(b) respectively. Reconstruction is obtained in 33 iterations. The center value of the inhomogeneities corresponding to Fig. 8(a) and 8(b) is 0.019 mm^{-1} and 0.13 mm^{-1} respectively. Horizontal cross-sections at y=41 through the centre of reconstructed images are presented in Fig. 9.

(a).Horizontal cross-sections at y=41 mm through the centre of images Fig. 3(b) and 6(a).

(b).Horizontal cross-sections at y=41 mm through the centre of images Fig. 4(b) and 7.

(c).Horizontal cross-sections at y=41 mm through the centre of images Fig. 5(c) and 8(a).

(d).Horizontal cross-sections at y=41 mm through the centre of images Fig. 5(d) and 8(b).

Finally to bring out the effectiveness of pMC-based algorithm, experimental data is used in a diffusion equation-based algorithm, which did not converge to meaningful results. Even with the 2% noise added simulation data, the DE-based algorithm has failed to converge. Without any noise in the simulated integrated intensity, as a test case, the object with two inhomogeneities (Fig. 5(a) and 5(b)) are reconstructed and the results are presented in Fig. 10(a) and 10(b). Even here the reconstruction was carried out with the same a priori information as in Fig. 5(c)and 5(d). The centre of the reconstructed inhomogeneities values in µ_{a}- (Fig. 10(a)) and µ_{s}-distribution (Fig. 10(b)) are 0.011 mm^{-1} and 0.06 mm^{-1} respectively. Reconstructed µs-inhomogeneity (Fig. 10(b)) is falsely shown in the reconstruction result (original is as shown in Fig. 5(b)).

## 5. Discussion and conclusions

With these experimental results a confirmation of the usefulness of a Monte-Carlo-based method to reconstruct low-scattering objects, where the diffusion equation fails to hold, is done. The novelty of this method lies in getting the derivative information using a simple analytical equation (Eq. 1) and using the same in subsequent iterations. Further, experimentally, it is shown that cross talk between absorption- and scattering coefficients is negligible in the pMC-based reconstructions. The DE-based reconstructions are prone to cross-talk as the diffusion- and absorption coefficients are not independent of each other. Therefore µ_{a}- and µ_{s} reconstructions based on sensitivity matrices calculated with respect to absorption and diffusion coefficients will have cross-talk which needs to be minimized by adding correction terms to the sensitivity matrix [33,34]. In addition, it is also seen that with the DE a measurement data type like intensity has a preferential weightage towards reconstructing µ_{a} and has poor sensitivity towards µ_{s} changes [34]. This is not so with an RTE or MC-based method as demonstrated through simulations and experiments here. In comparison, a large false positive is seen in the µ_{s} image at the µ_{a}-inhomogeneity location in the DE-based reconstruction shown in Fig. 10. Usefulness of dividing the reconstruction into a number of sub-problems based on angle of projection (view) in reducing computational complexity is demonstrated which give accurate reconstructions when location of inhomogeneities is *a priori* known.

The main applications of optical tomography being brain and breast imaging [1–7], a great challenge in the medical optical imaging is to develop a model which can reconstruct optical properties with a large range of variation. An example is the low-scattering CSF layer surrounding a highly scattering inhomogeneous region. To generalize, the typical problem would be to reconstruct low-scattering inclusions in an otherwise high scattering background. Whereas the RTE, and its equivalent the MC simulation, would hold good for both high and low-scattering regions, the DE is valid only where scattering predominates. One strategy would be to use the DE, which is computationally efficient to implement and switch to the RTE when one encounter the inclusions. This can be successful only if one can provide the correct boundary and source conditions at the interface between the inclusions and the background, which in itself is complicated. Another alternative approach can be to use MC as the forward model and pMC to provide update to both the Jacobian and the computed forward data which is computationally feasible. Reconstruction is limited only to the regions of interest, namely, the low-scattering inclusions, which are *a priori* obtained by analyzing the intensity data gathered around the objects following the procedure suggested in Ref. [43]. Identification of the regions of interest will drastically reduce the problem dimension, application of pMC to update the sensitivity matrix and the forward data will render computational cost manageable and the MC simulation which takes into account angles of scattering with its larger applicability will reconstruct low-scattering inclusions better, i.e., more accurately and with lesser inter-parameter cross-talk. The future direction for this work will be to demonstrate, through simulations and experiment, imaging of low-scattering inclusions in a high scattering background using the above strategy, which would justify the practical utility of the method presented here.

In Ref. [44] a three-dimensional (3D) Monte-Carlo code has been presented for modeling photon transport in adult head. Availability of such efficient 3D models pave the way for extending the 2D demonstration presented here to 3D imaging and using it for practical brain tissue characterization. This method could be combined with and guided by other imaging modalities such as MRI, X-ray CT and ultrasound to complement its capabilities to realize a superior imaging modality which gives more accurate and useful spatial as well as structural information of the human brain than what each one of them can individually provide.

## Acknowledgments

Authors are grateful to Prof. Hamid Dehghani for their useful discussions on diffusion-based methods. P.K.Y. is grateful to the Chairman, Supercomputer Education and Research Centre (SERC), Indian Institute of Science (IISc), Bangalore for providing the computational facilities for carrying out this work. P.K.Y.’s current address: Thayer School of Engineering, Dartmouth College, 8000 Cummings Hall, Hanover, NH 03755, USA. e-mail: Phaneendra.K.Yalavarthy@dartmouth.edu

## References and links

**1. **J. C. Hebden, D. A. Boas, J. S. George, and Anthony J. Durkin, “Topics in Biomedical Optics: Introduction,” Appl. Opt. **42**, 2869–2870 (2003). [CrossRef]

**2. **B. Chance, R. R. Alfano, B. J. Tromberg, M. Tamura, and E. M. Sevick-Muraca, ed., *Optical Tomography and Spectroscopy of Tissue IV*, Proc. SPIE4250 (2001).

**3. **B. Chance and R. R. Alfano, ed., *Optical tomography and spectroscopy of tissue: Theory, instrumentation, model, and human studies*, Proc. SPIE2979 (1997).

**4. **B. Chance and R. R. Alfano, ed., *Optical tomography, photon migration and spectroscopy of tissue and model media: theory, human studies, and instrumentation*, Proc. SPIE2389 (1995).

**5. **A. Yodh and B. Chance, “Spectroscopy and Imaging with diffusing light,” Phy. Today **48**, 34–40 (1995). [CrossRef]

**6. **J. C. Hebden, S. R. Arridge, and D. T. Delpy, “Optical imaging in medicine. I. Experimental techniques,” Phys. Med. Biol. **42**, 825–840 (1997). [CrossRef] [PubMed]

**7. **S. R. Arridge and J. C. Hebden, “Optical imaging in medicine. II. Modelling and reconstruction,” Phys. Med. Biol. **42**, 841–853 (1997). [CrossRef] [PubMed]

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

**9. **R. Chandrasekhar, *Radiation Transfer* (Oxford, Clarendon, 1950).

**10. **A. Ishimaru, *Wave Propagation and Scattering in Random Media* (IEEE Press, New York, 1997).

**11. **M. Firbank, S. R. Arridge, M. Schweiger, and D. T. Delpy, “An investigation of light transport through scattering bodies with nonscattering regions,” Phys. Med. Biol. **41**, 767–783 (1996). [CrossRef] [PubMed]

**12. **J. Ripoll, N. Nieto-Vesperinas, S. R. Arridge, and H. Dehghani, “Boundary conditions for light propagation in diffuse media with nonscattering regions,” J. Opt. Soc. Am. A **17**, 1671–1681 (2000). [CrossRef]

**13. **H. Dehghani, S. R. Arridge, M. Schweiger, and D. T. Delpy, “Optical tomography in the presence of void regions,” J. Opt. Soc. Am. A **17**, 1659–1670 (2000). [CrossRef]

**14. **Y. Xu, Q. Zhang, and H. Jiang, “Optical image reconstruction of non-scattering and low scattering heterogeneities in turbid media based on the diffusion approximation model,” J. Opt. A: Pure Appl. Opt. **6**, 29–35 (2004). [CrossRef]

**15. **E. Okada, M. Firbank, M. Schweiger, S. R. Arridge, M. Cope, and D. T. Delpy, “Theoretical and experimental investigation of near-infrared light propagation in a model of the adult head,” Appl. Opt. **36**, 21–31 (1997). [CrossRef] [PubMed]

**16. **M. Wolf, M. Keel, V. Dietz, K. von Siebenthal, H. U. Bucher, and O. Baenziger, “The influence of a clear layer on near-infrared spectrophotometry measurements using a liquid neonatal head phantom,” Phys. Med. Biol. **44**, 1743–1753 (1999). [CrossRef] [PubMed]

**17. **H. Dehghani and D. T. Delpy, “Near-infrared spectroscopy of the adult head: effect of scattering and absorbing obstructions in the cerebrospinal fluid layer on light distribution in the tissue,” Appl. Opt. **39**, 4721–4729 (2000). [CrossRef]

**18. **E. Okada and D. T. Delpy, “Near-infrared light propagation in an adult head model. I. Modeling of low-level scattering in the cerebrospinal fluid layer,” Appl. Opt. **42**, 2906–2914 (2003). [CrossRef] [PubMed]

**19. **E. Okada and D. T. Delpy, “Near-infrared light propagation in an adult head model. II. Effect of superficial tissue thickness on the sensitivity of the near-infrared spectroscopy signal,” Appl. Opt. **42**, 2915–2922 (2003). [CrossRef] [PubMed]

**20. **A. D. Klose, V. Prapavat, O. Minet, J. Beuthan, and G. Muller, “RA diagnostics applying optical tomography in frequency-domain,” Proc. SPIE **3196**, 194–204 (1997). [CrossRef]

**21. **A. D. Klose, A. H. Hielscher, K. M. Hanson, and J. Beuthan, “Three-dimensional optical tomography of a finger joint model for diagnostic of rheumatoid arthritis,” Proc. SPIE **3566**, 151–159(1998). [CrossRef]

**22. **Y. Xu, N. V. Iftimia, H. Jiang, L. L. Key, and M. B. Bolster, “Imaging of in vitro and in vivo bones and joints with continuous-wave diffuse optical tomography,” Opt. Express **8**, 447–451 (2001), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-8-7-447 [CrossRef] [PubMed]

**23. **O. Dorn, “A transport-backtransport method for optical tomography,” Inverse Problems **14**, 1107–1130 (1998). [CrossRef]

**24. **A. D. Klose and A. H. Hielscher, “Iterative reconstruction scheme for optical tomography based on the equation of radiative transfer,” Med. Phys. **26**, 1698–1707 (1999). [CrossRef] [PubMed]

**25. **A. D. Klose, U. Netz, J. Beuthan, and A. H. Hielscher, “Optical tomography using the time-independent equation of radiative transfer - Part 1:forward model,” J. Quant. Spectrosc. Radiat. Transf. **72**, 691–713 (2002). [CrossRef]

**26. **A. D. Klose and A. H. Hielscher, “Optical tomography using the time-independent equation of radiative transfer - Part 2:inverse model,” J. Quant. Spectrosc. Radiat. Transf. **72**, 715–732 (2002). [CrossRef]

**27. **A. D. Klose and A. H. Hielscher, “Quasi-Newton methods in optical tomographic image reconstruction,” Inverse Problems **14**, 387–403 (2003). [CrossRef]

**28. **M. Firbank, E. Okada, and D. T. Delpy, “A theoretical study of the signal contribution of regions of the adult head to near infrared spectroscopy studies of visual evoked responses,” Neuroimage **8**, 69–78 (1998). [CrossRef] [PubMed]

**29. **S. R. Arridge, H. Dehghani, M. Schweiger, and E. Okada, “The finite element model for the propagation of light scattering media: a direct method for domain with nonscattering regions,” Med. Phys. **27**, 252–264 (2000). [CrossRef] [PubMed]

**30. **H. Jiang, “Optical image reconstruction based on the third-order diffusion equations,” Opt. Express **4**, 241–246 (1999), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-4-8-241 [CrossRef] [PubMed]

**31. **N. Herschkowitz, “Brain development in the fetus, neonate and infant, Biol. Neonate **54**, 1–19 (1988) [CrossRef] [PubMed]

**32. **S. R. Arridge and W. R. B. Lionheart, “Nonuniqueness in diffusion-based optical tomography,” Opt. Lett. **23**, 882–884 (1998). [CrossRef]

**33. **Y. Pei, H. L. Graber, and R. L. Barbour, “Normalized-constraint algorithm for minimizing inter-parameter crosstalk in DC optical tomography,” Opt. Express **9**, 97–109 (2001), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-9-2-97 [CrossRef] [PubMed]

**34. **J. C. Ye, K. J. Webb, R. P. Millane, and T. J. Downar, “Modified distorted Born iterative method with an approximate Frechet derivative for optical diffusion tomography,” J. Opt. Soc. Am. A **16**, 1814–1826 (1999). [CrossRef]

**35. **A. Sassaroli, C. Blumetti, F. Martelli, L. Alianelli, D. Contini, A. Ismaelli, and G. Zaccanti, “Monte Carlo procedure for investigating light propagation and imaging of highly scattering media,” Appl. Opt. **37**, 7392–7400 (1998). [CrossRef]

**36. **C. K. Hayakawa, J. Spanier, F. Bevilacqua, A. K. Dunn, J. S. You, B. J. Tromberg, and V. Venugopalan, “Perturbation Monte Carlo methods to solve inverse photon migration problems in heterogeneous tissues,” Opt. Lett. **26**, 1335–1337 (2001). [CrossRef]

**37. **Y. Phaneendra Kumar and R. M. Vasu, “Reconstruction of optical properties of low-scattering tissue using derivative estimated through perturbation Monte-Carlo method,” J. Biomed. Opt. **9**, 1002–1012 (2004). [CrossRef] [PubMed]

**38. **F. Natterer and F. Wübbeling, “A propagation-backpropagation method for ultrasound tomography,” Inverse Problems **11**, 1225–1232 (1995). [CrossRef]

**39. **S. R. Arridge, M. Schweiger, M. Hiraoka, and D. T. Delpy, “A finite element approach for modeling photon transport in tissue,” Med. Phys. **20**, 299–309 (1993). [CrossRef] [PubMed]

**40. **A. H. Hielscher, A. D. Klose, and K. M. Hanson, “Gradient-based iterative image reconstruction scheme for time-resolved optical tomography,” IEEE Trans. Med. Imag. **18**, 262–271 (1999). [CrossRef]

**41. **J. Spanier and E. M. Gelbard, *Monte Carlo Principles and Neutron Transport Problems* (Addison-Wesley, Reading, Mass., 1969).

**42. **B. Jain, P. K. Gupta, V. A. Podzyavnikov, and V. K. Chevokin, “Development and characterization of a UV-visible streak camera and its use for time resolved fluorescence studies on human tissues,” Proc. National Laser Symposium, B.A.R.C., Mumbai, India, January 17–19, 1996, National Laser Program, Department of Atomic Energy, Government of India, C3–C4 (1996).

**43. **B. Kanmani and R. M. Vasu, “Diffuse optical tomography using intensity measurements and the *a priori* acquired regions of interest: theory and simulations,” Phy. Med. Biol. **50**, 247–264 (2005). [CrossRef]

**44. **D. A. Boas, J. P. Culver, J. J. Stott, and A. K. Dunn, “Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head,” Opt. Express **10**, 159–170 (2002), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-3-159 [PubMed]