## Abstract

We describe a novel Monte Carlo code for photon migration through 3D media with spatially varying optical properties. The code is validated against analytic solutions of the photon diffusion equation for semi-infinite homogeneous media. The code is also cross-validated for photon migration through a slab with an absorbing heterogeneity. A demonstration of the utility of the code is provided by showing time-resolved photon migration through a human head. This code, known as ‘tMCimg’, is available on the web and can serve as a resource for solving the forward problem for complex 3D structural data obtained by MRI or CT.

© 2002 Optical Society of America

## 1. Introduction

Imaging tissue with diffuse light (i.e. diffuse optical tomography) produces images with poor structural details (for instance non-invasive measurements of the human cortex have a resolution typically no better than 5–10 mm), but the images contain exquisite functional information that complements the information obtained by other imaging modalities such as MRI, X-Ray CT, and PET [1]. Because of limited spatial resolution, diffuse optical tomography is increasingly being used in combination with imaging methods that provide high-resolution structural information, such as MRI, CT, X-Ray, and ultrasound [2–4]. To incorporate high-resolution structural information into the diffuse optical tomography problem it is necessary to have accurate forward models for highly heterogeneous medium with arbitrary boundaries. We present a computer code for accurate forward solutions of complex 3D media and demonstrate its utility using a 3D MRI segmentation of the human head as the structure for photon migration.

During the last 15 years, the modeling of photon migration within tissue has allowed for the quantitative determination of tissue hemoglobin oxygen saturation for basic physiological research and clinical applications [1, 5, 6]. This progress has motivated the extension of photon migration techniques to tissue imaging applications using diffuse optical tomography (DOT) methods. These imaging methods enable determination of spatially varying functional and structural information. More accurate mapping of tissue will increase the utility of photon migration technology for a broader range of applications including: optical breast imaging [7–9], functional brain imaging [6, 10–12], as well as imaging stroke [6] and head trauma [13, 14]. Progress in diffuse optical tomography research has been promising. Spatial variations in absorption and scattering coefficients can be experimentally imaged under a wide variety of conditions in simplified tissue phantoms with piece-wise continuous optical properties [15, 16]. In addition, some groups have demonstrated the possibility of imaging within more realistic head models [2, 17, 18].

Improving the accuracy of diffuse optical tomography of complex tissue structures depends on the development of more sophisticated methods for solving the forward photon migration problem. While finite-element and finite-difference solutions of the diffusion equation are likely to provide accurate solutions, they need to be compared against solutions of the radiative transport equation [19, 20] to verify accuracy near internal and external boundaries. Furthermore, there are relevant conditions under which the diffusion approximation introduces significant errors including highly absorbing regions (e.g. brain bleeds) and weakly scattering regions (such as cerebral spinal fluid in the brain) [21–23]. Hielscher et al [24] have implemented a finite-difference solution of the radiative transport equation and compared it against solutions of the diffusion equation for homogeneous and complex heterogeneous media like the human head. Arridge et al [25, 26] have developed a hybrid diffusion - radiosity method for handling nonscattering voids within highly scattering media. Dorn has developed an efficient approach for utilizing the transport equation in the tomography problem[27]. While these approaches hold great promise, they are complex and there are several issues involved with establishing their accuracy (e.g. the number of nodes used to discretize the boundaries, spatial and angular coordinates). Monte Carlo solution methods, on the other hand, are conceptually simpler to implement and rely on fewer assumptions, but at the expense of computational speed.

We describe a novel Monte Carlo code that allows us to obtain results, for example, in a complex 3D head model with a signal-to-noise ratio greater than 100 up to distances of 30 mm with a 1 mm^{2} detector with 10^{8} photons propagated within 5–10 hours of computer time on a Pentium III 1000 MHz CPU. The method is validated against an accepted analytic solution for a semi-infinite medium, and cross-validated for a slab geometry with an absorbing inclusion. We then illustrate the utility of the Monte Carlo method with a novel simulation and movie of time-resolved photon migration through a human head and deduce the depth sensitivity of different measurement types. The human head model was obtained from a 3D segmented anatomical MRI. The ability to perform such simulations in a medium with arbitrary boundaries and spatial variation in the optical properties, in our publicly available code ‘tMCimg’ [28], is a significant advancement over the capabilities of the widely used and publicly available code known as ‘MCML’ developed by Wang and Jacques [29] which models photon propagation through layered media with planar boundaries.

## 2. Method

The rules used in Monte Carlo programs for photon migration in biological tissue are straightforward as described by Wang and Jacques [29, 30]. Differences between programs arise in how the photon fluence or flux is recorded during the Monte Carlo simulation. Below, we briefly describe the rules for photon migration and how we record the photon fluence within the medium and the flux of photons exiting the medium.

To begin, the initial position and direction of the photon is defined. In our case, we consider a point source of photons with a defined direction, entering the medium on the surface. The source can be diverging, as defined by the source numerical aperture (NA), in which case the azimuthal angle is determined by a random number uniformly distributed from 0 to 2π and the elevation (angle from the source axis) is determined from another random number distributed uniformly from 0 to sin^{-1}(NA).

Given the initial position and direction of the photon, the length to the first scattering event is calculated from an exponential distribution. Absorption of the photon is considered by decreasing the photon weight by exp(-μ_{a}L) where μ_{a} is the absorption coefficient and L is the length traveled by the photon [29]. The photon is moved this length. A scattering angle is calculated using the probability distribution given by the Henyey-Greenstein phase function [29]. A new scattering length is then determined from an exponential distribution. The photon is propagated this new length in the new direction. This process continues until the photon exits the medium or has traveled longer than a predefined period of time. We typically stop the photon after 10 ns since the probability of photon detection in tissue after such a period of time is exceedingly small. When the photon does try to leave the medium, the probability of an internal reflection is calculated using Fresnel’s equation [29, 31]. If a reflection occurs, then the photon is reflected back into the medium the appropriate distance and migration continues. Otherwise, the migration of that particular photon halts and a new photon is launched into the medium at the predefined source location.

We consider photon migration through media with spatially varying optical properties by representing the optical properties within voxels on a cubic grid. We typically consider a grid spacing of 1 mm, but this can be scaled. As the photon is propagated from one scattering event to the next, a check is made every 1 grid spacing to see if the scattering or absorption coefficient has changed. If the scattering coefficient has changed, the remainder of the scattering length is renormalized by ${{\mu}_{s}}^{\mathit{\text{old}}}$
/${{\mu}_{s}}^{\mathit{\text{new}}}$
where *μ*_{s}
is old and new scattering coefficient, as described in Jacques and Wang [32]. If the absorption coefficient has changed, then the photon weight is decreased with the new absorption coefficient [32]. The scattering angle is determined by the value of the scattering anisotropy factor, *g*, within the particular voxel containing the scattering event.

The photon fluence within the medium is determined by accumulating the photon weight every 1 grid spacing in the voxel corresponding to the present position of the photon. When a photon exits the scattering medium and enters the surrounding non-scattering medium, the exiting photon flux is determined by accumulating the photon weight in a bin representing the non-scattering voxel that the photon entered. In addition, if the photon exits in a location that has been predefined as a detector location, then information is recorded in a history file identifying the detector and pathlength of the photon in each tissue type prior to being detected. The tissue types are specified as input to the Monte Carlo code and refer to voxel groups with identical optical properties. The code supports up to 100 different tissue types. This information can be used in post-processing to analyze the effect of absorption changes within different tissue types as described further below.

After all simulated photons have been propagated, it is necessary to normalize the exiting photon flux and the photon fluence within the medium. The exiting photon flux , *J*_{out}
(**r**), is divided by the number of simulated photons. Normalizing the photon fluence, *Φ*(**r**), is more involved. To conserve energy, the exiting photon flux plus the number of photons absorbed in the medium must equal the number of simulated photons, which we normalize to 1. The number of photons absorbed at a given point is *Φ*(**r**) μ_{a}(**r**). Therefore,

where **r**
_{i} indicates the position of each surface and volume voxel, *V*_{voxel}
is the voxel volume, and *A*_{i}
is the area of the surface element at position **r**
_{i}. The normalization factor for *Φ*(**r**) is determined from this relation.

We now described how the detected photons stored in the history file are analyzed. First, it is necessary to realize that the photon fluence can be written in the form [33]

where Φ(*t*) is the measured photon fluence, *N*_{photons}
(*t*) is the number of photons collected in a time-gate of width ∆*t* centered at time *t*, exp(-*μ*_{aj}*L*_{ij}
) accounts for the effects of absorption in each region where *L*_{ij}
is the pathlength of photon *i* through region *j*, and the photon migration time is related to the photon pathlength by the speed of light in the medium. The history file of detected photons provides all the necessary information for calculation of eq. (2).

## 3. Solutions of the Diffusion Equation for Comparison with Monte Carlo

The solution of the photon diffusion equation [19, 31, 34] for a semi-infinite homogeneous medium for a continuous-wave point source of light is given by [35, 36]

and for a time-domain measurement with an impulse point source of light [34, 36]

Φ(**r**
_{s},**r**
_{d},*t*) is the photon fluence at the detector position **r**
_{d} at time t, generated by a point source of amplitude *S* at position **r**
_{s}. *D*=*v*/(*3μ*_{s}*’*) is the photon diffusion coefficient [37, 38], *μ*_{s}*’* is the reduced scattering coefficient, *μ*_{a}
is the absorption coefficient, and *v* is the speed of light in the medium. The semi-infinite boundary condition is satisfied by the method of images [31, 36]. The position of the image source is indicated by **r**
_{s,i}.

When the optical properties are spatially varying, there are two standard approaches to finding approximate linear solutions: the first Born approximation [19, 39]

and the Rytov approximation [39]

Each approach decomposes the total fluence Φ into Φ_{o} which only depends on the background optical properties *μ*_{ao}
and *μ*_{so}*’*, and Φ_{pert} which is linearly related to spatial variations in the optical properties *δμ*_{a}
and *δμ*_{s}*’*. Historically, the Born approximation is more common, but for the diffuse optical tomography problem the Rytov approximation tends to be more accurate as it accounts for some non-linear saturation due to increasing perturbation in the absorption coefficient [10]. For the Rytov approximation, assuming absorption variations only,

**r**
_{s} and **r**
_{d} are the position of the source and detector respectively, *G* is the Greens function of the photon diffusion equation for the background optical properties given the boundary conditions. For the Born approximation, Φ_{pert} is not normalized by Φ_{o}(**r**
_{s},**r**
_{d}) If the background is homogeneous then Φ_{o} and *G* can be expressed analytically in some simple geometries [34, 40], where the difference between Φ_{o} and *G* is amplitude factor *S* (see eq. (4)).

#### Validation of the Monte Carlo Code in a Homogeneous Semi-Infinite Medium

The first step in validating the Monte Carlo code is to compare the results for the remitted flux of photons from a semi-infinite homogeneous medium with an analytic solution of the diffusion equation as has been done in the past by several others [29, 36]. For the simulation, a collimated point source was incident on a semi-infinite medium with a scattering coefficient of 1 mm^{-1}, a scattering anisotropy of *g* = 0.01, and an absorption coefficient of 0.005 mm^{-1}. The medium had dimensions of 60 × 60 × 60 mm with the source positioned at (x,y,z) = (30,30,1) mm. All boundaries were treated as index matched. The source was sufficiently far from the edges so that the medium is effectively semi-infinite. A simulation was executed with 10^{8} photons which took approximately 6 hours on a 1 GHz Pentium 3.

The results for the continuous-wave simulation are shown in fig. 1. Fig. 1a compares the radially resolved flux of photons exiting the medium determined by Monte Carlo (eq. (2)) and eq. (3), while fig. 1b compares the fluence field within the medium with contours drawn every 10 dB. In both cases we observe decent agreement. In fig. 1a, the systematically smaller Monte Carlo result at larger separations results from the influence of the edge at 30 mm which is not modeled in the solution of the diffusion equation. The discrepancy observed in the fluence field in fig. 1b results from the isotropic point source approximation, within the diffusion equation, of the collimated source. This isotropic point source approximation gives rise to decent agreement for the remitted flux of photons on the surface (as shown in fig. 1a), but at the expense of decent agreement within the highly scattering medium, especially in the forward direction. Better agreement is found within the medium by extending the isotropic point source slightly deeper into the medium, but at the expense of reduced agreement for the remitted flux. A more accurate solution is to model the exponentially attenuated collimated source within the diffuse equation [31].

The results for the time-domain simulation are shown in fig. 2. In fig. 2a the time-response profile is shown for the remitted flux at x = 15 mm from the source and for the photon fluence at position (x, z) = (15, 10) mm within the medium, for both the Monte Carlo simulation and the diffusion equation. Fig. 2b compares the spatial distribution of the photon fluence within the medium by plotting the half-maximum contours at time points of 0.1 ns, 0.5, 1.0, 1.5, and 2.0 ns. The temporal and spatial agreement is fairly good. The systematic decrease for the Monte Carlo photon fluence within the medium at longer times results from the edge effects. At 0.1 ns the Monte Carlo spatial distribution is much more compact than that for the diffusion equation as expected from the non-casual behavior of the diffusion equation.

## 4. Validating the Monte Carlo code in a slab geometry with an absorbing inclusion

We have three approaches for analyzing the effect of an absorbing inclusion with different absorption coefficients on the detection of diffuse photons: 1) run a Monte Carlo simulation for each absorption coefficient, 2) utilize the history file of detected photons which records the pathlength each detected photon spent in each piece-wise constant region, and 3) utilize the Born and/or Rytov approximation which can be produced from the product of the fluence field from a source into the medium, *Φ*_{o}
(**r**
_{s},**r**), and from the detector into the medium, *G*(**r**,**r**
_{d}), (see eq. (7), so called adjoint fields [41, 42]). Comparing these three approaches serves to cross-validate the methodology for media with spatially varying optical properties.

The comparison is performed in a 40 mm thick slab with a cubic absorbing inclusion 15 × 15 × 15 mm centered 17 mm into the slab from the source. The source and detector are positioned directly across each other centered on the inclusion (see fig. 3a). The scattering coefficient of the medium was set to 1 mm^{-1}, the scattering anisotropy factor was set to 0.01, the absorption coefficient of the background was set to 0.005 mm^{-1}, while that of the inclusion was increased from 0.005 mm^{-1} to 0.065 mm^{-1}. A full Monte Carlo simulation (which traces out the scattering path of every photon) was executed with the different absorption coefficients to calculate the flux of photons received at the detector. These results normalized by the flux detected with no absorbing inclusion are shown by the square symbols in fig. 3b. Notice that as the absorption coefficient initially increases, there is an approximate linear decrease in the intensity which then saturates for higher absorption coefficients. The detected flux of photons for different absorption coefficients can also be calculated from the photon history file produced by a single Monte Carlo simulation, as described above in eq. (2).

These results are shown by the solid line in fig. 3b. As expected, the results produced by the history file (generated by a previous full Monte Carlo simulation but using a different value of the absorption coefficient in the post-processing) agree with full Monte Carlo results. Finally, the results produced by the Born and Rytov approximations are shown by the dashed and dot-dashed lines respectively. These approximations are linear and exponential respectively and therefore show the expected agreement for small changes in the absorption coefficient, but the approximations break down for larger absorption coefficients as evidenced by the poor agreement with the non-linear Monte Carlo results. Notice that the Rytov approximation shows the expected better agreement with the Monte Carlo results. Finally, in fig. 4 we show the photon fluence within the medium without the absorber and the change with the absorbing inclusion demonstrating the shadowing effect of the absorbing inclusion.

## 5. Full 3D Head

Figure 5 shows an anatomical MRI of a human head, segmented into five tissue types (air, scalp, skull, cerebral spinal fluid, and gray/white matter; see [43]), with a contour overlay indicating the photon migration spatial sensitivity profile for (a) continuous-wave, (b) 200 MHz modulation, and (c,d) pulsed measurements. The head was contained within 151×171×232 1mm^{3} voxels and the Monte Carlo simulation recorded the temporal response and took approximately 10 hours. The continuous-wave and 200 MHz results were obtained from the modulus of Fourier transform of the temporal response. One contour line is shown for each half order of magnitude (10 dB) signal loss, and the contours end after 3 orders of magnitude in loss (60 dB). For the 3D Monte Carlo simulation, we assumed that μ_{s}’ = 1 mm^{-1} and μ_{a}=0.04 mm^{-1} for the scalp and skull, μ_{s}’ = 0.01 mm^{-1} and μ_{a}=0.001 mm^{-1} for the CSF, and μ_{s}’ = 1.25 mm^{-1} and μ_{a}=0.025 mm^{-1} for the gray/white matter. Note how the contours extend several millimeters into the brain tissue, indicating sensitivity to changes in cortical optical properties. The depth penetration difference between the continuous-wave and 200 MHz measurements is difficult to discern. A ratio of the two sensitivity profiles (not shown) shows that the 200 MHz profile is shifted slightly towards the surface. The time-domain sensitivity profiles suggest the possibility of obtaining greater penetration depths in the head from measurements made at longer delay times. This is further supported by the movie of the temporal evolution of a light pulse within the head as shown in fig. 6. This movie illustrates the usefulness of the code for visualizing the temporal evolution of photon migration through a heterogeneous medium.

In fig. 7 we show an estimate of the signal-to-noise ratio obtained by running 10^{8} photons on the 3D human head, as determined from running 9 independent Monte Carlo simulations with different random number seeds. Fig. 7a shows the flux of photons exiting from the head as a function of distance from the source voxel to each voxel into which photons escaped, i.e. voxels describing the air surrounding the head. This flux is normalized by the number of simulated photons and of course depends on the 3D geometry and spatially varying optical properties. The expected exponential decay of the photon flux with distance is observed. The structure is seen in the data because of the cross-sectional area of different air voxels against the head. Some air voxels adjoined the head on 1, 2, or more surfaces, while some voxels only touched at the corner. The needed correction factor is just the effective cross-sectional area of each air voxel (see eq. (1)). At present, the code does not correct for this cross-sectional area for a curved surface. The noise was determined from the standard deviation of the 9 independent Monte Carlo simulations. This deviation is seen in fig. 7a by the different color circles not perfectly overlapping at larger separations. The signal-to-noise ratio versus separation is shown in fig. 7b.

Note that an SNR greater than 100 is achieved for separations smaller than 2 cm, while an SNR greater than 10 is still found for separations up to 5 cm. The SNR can be improved by the square-root of the area of the detector by averaging over neighboring voxels. The voxel size in these measurements were 1 mm. Typical detector sizes for human head measurements are 3–5 mm. Using detectors of this size will increase the SNR by a factor of 3–5. The SNR will also improve as the square-root of the number of simulated photons. Finally, the SNR will also improve by decreasing the reduced scattering coefficient and/or the absorption coefficient so that a greater number of photons are able to travel greater distances through the tissue (note that reducing the scattering coefficient can actually decreasing the SNR within a few millimeters of the source).

## 6. Summary

We have described a computer code for calculating the migration of light through 3D highly scattering media with spatially varying optical properties and arbitrary boundary conditions. As diffuse optical tomography is combined with and guided by other imaging modalities with superior spatial resolution (and thus superior structural information), it becomes increasingly important to have an accurate photon migration forward solver. This Monte Carlo computer code can take structural information provided by MRI, or X-Ray CT for example, and accurately solve the photon migration forward problem in a reasonable amount of time. In addition to allowing validation of more efficient forward solvers, this code should prove useful on its own. As a few examples, this code could be used for investigating the feasibility of quantitative tissue characterization with different data types, cerebral oximetry with time domain data for example, aiding in the development of an optimal source-detector geometry for imaging brain function, as well as serving as the forward solver for MRI guided diffuse optical tomography of the head.

DAB acknowledges funding from Advanced Research Technologies, NIH R29-NS38842, NIH P41-RR14075 and from the Center for Innovative Minimally Invasive Therapies. This research was funded in part by the US Army, under Cooperative Agreement No. DAMD17-99-2-9001. The material presented does not necessarily reflect the position or the policy of the Government, and no official endorsement should be inferred. A.K.D. acknowledges funding from NIH K25 NS41291-01.

## References and Links

**1. **A. Villringer and B. Chance, “Non-invasive optical spectroscopy and imaging of human brain function,” Trends Neurosci. **20**, 435–442 (1997). [CrossRef] [PubMed]

**2. **B. W. Pogue and K. D. Paulsen, “High-resolution near-infrared tomographic imaging simulations of the rat cranium by use of a priori magnetic resonance imaging structural information,” Opt. Lett. **23**, 1716–1718 (1998). [CrossRef]

**3. **Q. Zhu, T. Durduran, V. Ntziachristos, M. Holboke, and A. G. Yodh, “Imager that combines near-infrared diffusive light and ultrasound,” Opt. Lett. **24**, 1050–1052 (1999). [CrossRef]

**4. **M. J. Holboke, B. J. Tromberg, X. Li, N. Shah, J. Fishkin, Kidney D., J. Butler, B. Chance, and A. G. Yodh, “Three-dimensional diffuse optical mammography with ultrasound localization in a human subject,” J Biomed Opt **5**, 237–47. (2000). [CrossRef] [PubMed]

**5. **D. A. Benaron, W. F. Cheong, and D. K. Stevenson, “Tissue Optics,” Science **276**, 2002–2003 (1997). [CrossRef] [PubMed]

**6. **D. A. Benaron, S. R. Hintz, A. Villringer, D. Boas, A. Kleinschmidt, J. Frahm, C. Hirth, H. Obrig, J. C. van Houten, E. L. Kermit, W. F. Cheong, and D. K. Stevenson, “Noninvasive functional imaging of human brain using light,” J Cereb Blood Flow Metab **20**, 469–77 (2000). [CrossRef] [PubMed]

**7. **B. W. Pogue, S. P. Poplack, T. O. McBride, W. A. Wells, K. S. Osterman, U. L. Osterberg, and K. D. Paulsen, “Quantitative hemoglobin tomography with diffuse near-infrared spectroscopy: pilot results in the breast,” Radiology **218**, 261–6. (2001). [PubMed]

**8. **V. Ntziachristos and B. Chance, “Probing physiology and molecular function using optical imaging: applications to breast cancer,” Breast Cancer Res **3**, 41–6 (2001). [CrossRef] [PubMed]

**9. **V. Ntziachristos, A. G. Yodh, M. Schnall, and B. Chance, “Concurrent MRI and diffuse optical tomography of breast after indocyanine green enhancement,” Proc Natl Acad Sci U S A **97**, 2767–72 (2000). [CrossRef] [PubMed]

**10. **M. A. Franceschini, V. Toronov, M. Filiaci, E. Gratton, and S. Fanini, “On-line optical imaging of the human brain with 160-ms temporal resolution,” Opt. Express **6**, 49–57 (2000). http://www.opticsexpress.org/oearchive/source/18957.htm [CrossRef] [PubMed]

**11. **S. R. Hintz, D. A. Benaron, A. M. Siegel, A. Zourabian, D. K. Stevenson, and D. A. Boas, “Bedside functional imaging of the premature infant brain during passive motor activation,” J Perinat Med **29**, 335–43 (2001). [CrossRef] [PubMed]

**12. **A. Bluestone, G. Abdoulaev, C. Schmitz, R. Barbour, and A. Hielscher, “Three-dimensional optical tomography of hemodynamics in the human head,” Opt. Express **9**, 272–286 (2001). http://www.opticsexpress.org/oearchive/source/34858.htm [CrossRef] [PubMed]

**13. **C. S. Robertson, S. P. Gopinath, and B. Chance, “A new application for near-infrared spectroscopy: Detection of delayed intracranial hematomas after head injury,” Journal of Neurotrauma **12**, 591–600 (1995). [CrossRef] [PubMed]

**14. **S. R. Hintz, W. F. Cheong, J. P. van Houten, D. K. Stevenson, and D. A. Benaron, “Bedside imaging of intracranial hemorrhage in the neonate using light: comparison with ultrasound, computed tomography, and magnetic resonance imaging,” Pediatr Res **45**, 54–9. (1999). [CrossRef] [PubMed]

**15. **M. A. O’Leary, D. A. Boas, B. Chance, and A. G. Yodh, “Experimental images of heterogeneous turbid media by frequency-domain diffusing-photon tomography,” Opt. Lett. **20**, 426–428 (1995). [CrossRef]

**16. **J. C. Hebden, D. J. Hall, M. Firbank, and D. T. Delpy, “Time-resolved optical imaging of a solid tissue-equivalent phantom,“ Appl. Opt. **34**, 8038–8047 (1995). [CrossRef] [PubMed]

**17. **M. Schweiger and S. R. Arridge, “Optical tomographic reconstruction in a complex head model using a priori region boundary information,” Phys Med Biol **44**, 2703–21 (1999). [CrossRef] [PubMed]

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

**19. **A. Ishimaru, *Wave Propagation and Scattering in Random Media* (Academic Press, Inc., San Diego1978).

**20. **E. Okada, M. Schweiger, S. R. Arridge, M. Firbank, and D. T. Delpy, “Experimental validation of Monte Carlo and finite-element methods of estimation of the optical path length in inhomogeneous tissue,” Appl. Opt. **35**, 3362–3371 (1996). [CrossRef] [PubMed]

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

**22. **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]

**23. **M. Firbank, Okada E., 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]

**24. **A. H. Hielscher, R. E. Alcouffe, and R. L. Barbour, “Comparison of finite-difference transport and diffusion calculations for photon migration in homogeneous and heterogeneous tissues,” Phys Med Biol **43**, 1285–302. (1998). [CrossRef] [PubMed]

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

**26. **J. Ripoll, M. Nieto-Vesperinas, S. R. Arridge, and H. Dehghani, “Boundary conditions for light propagation in diffusive media with nonscattering regions,” J Opt Soc Am A Opt Image Sci Vis **17**, 1671–81. (2000). [CrossRef] [PubMed]

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

**28. **J. J. Stott and D. A. Boas, tMCimg: Monte Carlo code for photon migration through general 3D Media. http://www.nmr.mgh.harvard.edu/DOT

**29. **L. Wang, S. L. Jacques, and L. Zheng, “MCML-Monte Carlo modeling of light transport in multi-layered tissues,” Computer Methods and Programs in Biomedicine **47**, 131–146 (1995). [CrossRef] [PubMed]

**30. **S. L. Jacques and L. Wang, “Monte Carlo modeling of light transport in tissues” in *Optical-Thermal response of laser-irradiated tissue*, Welch and v. Gemert (Plenum, New York1995).

**31. **R. C. Haskell, L. O. Svaasand, T. Tsay, T. Feng, M. S. McAdams, and B. J. Tromberg, “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am A **11**, 2727–2741 (1994). [CrossRef]

**32. **S. L. Jacques and L. Wang, “Monte-Caro Modeling of Light Transport in Tissues” in *Optical-Thermal Response of Laser Irradiated Tissue*, A. J. Welch and M. C. J. van Gemert (Plenum, New York1995).

**33. **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]

**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**, 2331–2336 (1989). [CrossRef] [PubMed]

**35. **T. J. Farrell, M. S. Patterson, and B. Wilson, “A diffusion theory model of spatially resolved, steady state diffuse reflectance for the noninvasive determination of tissue optical properties in vivo,” Med. Phys. **19**, 879–888 (1992). [CrossRef] [PubMed]

**36. **A. Kienle and M. S. Patterson, “Improved solutions of the steady-state and the time-resolved diffusion equations for reflectance from semi-infinite turbid medium,” Journal of the Optical Society of America **14**, 246–254 (1997). [CrossRef] [PubMed]

**37. **K. Furutsu and Y. Yamada, “Diffusion approximation for a dissipative random medium and the applications,” Phys.Rev.E **50**, 3634 (1994). [CrossRef]

**38. **T. Durduran, B. Chance, A. G. Yodh, and D. A. Boas, “Does the photon diffusion coefficient depend on absorption?,” J. Opt. Soc. Am A **14**, 3358–3365 (1997). [CrossRef]

**39. **A. C. Kak and M. Slaney, *Principles of Computerized Tomographic Imaging* (IEEE Press, New York1988).

**40. **S. R. Arridge, M. Cope, and D. T. Delpy, “The theoretical basis for the determination of optical pathlengths in tissue: temporal and frequency analysis,” Phys Med Biol **37**, 1531–60 (1992). [CrossRef] [PubMed]

**41. **S. R. Arridge and M. Schweiger, “Photon-measurement density functions. Part2: Finite-element-method calculations,” Appl. Opt. **34**, 8026–8037 (1995). [CrossRef] [PubMed]

**42. **S. R. Arridge and M. Schweiger, “A gradient-based optimisation scheme for optical tomography,” Opt. Express **2**, 213–226 (1998). http://www.opticsexpress.org/oearchive/source/4014.htm [CrossRef] [PubMed]

**43. **A. M. Dale, B. Fischl, and M. I. Sereno, “Cortical surface-based analysis. I. Segmentation and surface reconstruction,” Neuroimage **9**, 179–94 (1999). [CrossRef] [PubMed]