Shape based Monte Carlo code for light transport in complex heterogeneous Tissues

doi:10.1364/OE.15.014086

« Show journal navigation

Shape based Monte Carlo code for light transport in complex heterogeneous Tissues

Eduardo Margallo-Balbás and Patrick J. French »View Author Affiliations

Optics Express, Vol. 15, Issue 21, pp. 14086-14098 (2007)
http://dx.doi.org/10.1364/OE.15.014086

View Full Text Article

Acrobat PDF (497 KB)

Journal Search
Article Lookup

Article Tools

Citations

Alert me when this article is cited
Export Citation/Save

Article Navigation

▸ Article Outline
– Abstract
1. Introduction
2. Implementation
3. Validation and performance
4. Application to trabecular bone
5. Conclusions
– References and links

Article Tools
Full Text
Article Info
References
Cited By
Figures
Metrics
Download PDF

Viewing Equations in the
HTML Article

Optimal Viewing Recommendations

Full Text
Article Info
References (45)
Cited By
Figures (6)
Metrics

Abstract

A Monte Carlo code for the calculation of light transport in heterogeneous scattering media is presented together with its validation. Triangle meshes are used to define the interfaces between different materials, in contrast with techniques based on individual volume elements. This approach allows to address realistic problems in a flexible way. A hierarchical spatial organisation enables a fast photon-surface intersection test. The application of the new environment to evaluate the impact of the trabecular structure of bone on its optical properties is demonstrated. A model of the trabecular micro structure recovered from microCT data was used to compute light distribution within tissue. Time-resolved curves across a spherical bone volume were computed. The work presented enables simulation of radiative transport in complex reality-based models of tissue and serves as a powerful, generic tool to study the effect of heterogeneity in the field of biomedical optics.

1. Introduction

In the last years the applications of optical methods in therapy and diagnosis have expanded significantly. A large amount of work has been performed on the utilisation of light to non-invasively study human tissues [1

X. Gu, Q. Zhang, M. Bartlett, L. Schutz, L. Fajardo, and H. Jiang, “Differentiation of cysts from solid tumours in the breast with Diffuse Optical Tomography,” Acad. Radiol. 11, 53–60 (2004). [CrossRef] [PubMed]

, 2

P. Taroni, G. Danesini, A. Torricelli, A. Pifferi, L. Spinelli, and R. Cubbedu, “Clinical trial of time-resolved scanning optical mammography at 4 wavelengths between 683 and 975 nm,” J. Biomed. Opt. 9(3), 464–473 (2004). [CrossRef] [PubMed]

] and to deliver various treatments [3

S. Brown, E. Brown, and I. Walker, “The present and future role of photodynamic therapy in cancer treatment,” Oncology 5, 497–508 (2004). [PubMed]

]. These positive results are based on a good understanding of light transport in turbid media, the availability of powerful numerical tools and a precise knowledge of the optical properties of the biological materials of interest.

Light transport in turbid media can be accurately modelled by means of the Radiative Transfer Equation (RTE) [4

S. Chandrasekhar, Radiative Transfer (Oxford University Press, London, 1960).

]. However, its integro-differential nature does not allow an analytical closed-form solution in practical situations and numerical methods and approximations are required to address realistic geometries. Monte Carlo (MC) methods represent a stochastic realisation of the RTE, having the potential to approximate the exact solution for a given problem with any desired accuracy. Additionally, they possess the ability to model realistic sources, detectors and sample geometries. Their main drawback is the need to simulate a large number of photon paths in order to reach a high accuracy level. Because of their benefits, many MC codes have been developed in the various fields of science and technology that deal with radiation transfer in turbid media. Besides biomedical optics [5

L. Wang, S. Jacques, and L. Zheng, “MCML - Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Meth. Prog. Bio. 47, 131–146 (1995). [CrossRef]

, 6

T. Pfefer, J. Barton, E. Chan, M. Ducros, B. Sorg, T. Milner, J. Nelson, and A. Welch, “A three-dimensional modular adaptable grid numerical model for light propagation during laser irradiation of skin tissue,” IEEE J. Sel. Top. Quantum Electron. 2(4), 934–942 (1996). [CrossRef]

, 7

D. Boas, J. Culver, J. Stott, and A. Dunn, “Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head,” Opt. Express 10(3), 159–170 (2002). [PubMed]

, 8

D. Côté and I. Vitkin, “Robust concentration determination of optically active molecules in turbid media with validated three-dimensional polarization sensitive Monte Carlo calculations,” Opt. Express 13, 148–163 (2005). [CrossRef] [PubMed]

], these fields include nuclear physics [9

F. Brown, R. Barrett, T. Booth, J. Bull, L. Cox, R. Forster, J. Goorley, R. Mosteller, S. Post, R. Prael, E. Selcow, A. Sood, and J. Sweezy, “MCNP Version 5,” Trans. Am. Nucl. Soc. 87, 273 (2002).

], radiative heat transfer [10

J. Howell, “The Monte Carlo method in radiative heat transfer,” J. Heat Transf. 120(3), 547–560 (1998). [CrossRef]

] and atmospheric optics [11

A. Marshak and A. B. Davis, eds., 3D Radiative Transfer in Cloudy Atmospheres (Springer-Verlag, Berlin, 2005). [CrossRef]

]. Each of the cited codes reflects the needs of the discipline for which it was developed in terms of problem modelling and magnitudes of interest.

In recent years, researchers in computer graphics have been working on physically realistic models of light transport for high quality rendering of virtual scenes [12

A. S. Glassner, Principles of Digital Image Synthesis (Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 1994).

, 13

E. Veach, “Robust Monte Carlo methods for light transport simulation,” Ph.D. thesis, Stanford University (1997).

]. Monte Carlo methods have been adapted to the needs of this field, sometimes including turbid media [14

E. Cerezo, F. Pérez, X. Pueyo, F. Serón, and F. Sillion, “A survey on participating media rendering techniques,” Visual. Comput. 21, 303–328 (2005). [CrossRef]

]. One of the differences between computer graphics and other research areas has been the focus on producing images [13

E. Veach, “Robust Monte Carlo methods for light transport simulation,” Ph.D. thesis, Stanford University (1997).

] rather than accurately computing some relevant physical magnitudes. Also the evaluation of a solution has different standards in this field: instead of trying to minimise the error relative to an exact solution based on an objective metric, these codes try to improve the scene quality as it is perceived by a human observer. Despite these differences, the emphasis on arbitrarily complex geometries, rather than idealised system models, has lead to important advances in efficient geometry organisation schemes [15

L. Szirmay-Kalos, V. Havran, B. Balász, and L. Szécsi, “On the efficiency of ray-shooting acceleration schemes,” in SCCG ’02: Proceedings of the 18th spring conference on Computer graphics , pp. 97–106 (ACM Press, New York, NY, USA, 2002). [CrossRef]

Since the work by Wang et al. [5

L. Wang, S. Jacques, and L. Zheng, “MCML - Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Meth. Prog. Bio. 47, 131–146 (1995). [CrossRef]

] for multi-layered tissues, there have been some important developments in MC codes able to deal with 3D geometries, especially the ones that can be recovered from imaging modalities like MRI or CT. The code described by Boas et al. is based on a voxelised model of the problem [7

] and it is similar to other works, like the one by Pfefer [6

]. Two problems are apparent in this type of volumetric description. First, the rigid space discretisation derived from fixed volume elements imposes a homogeneous resolution. This might not be appropriate, in terms of computational resources, if the problem shows heterogeneity at different dimensional scales. Second, a volume-based geometry description does not allow for explicit consideration of boundaries, which are important if the mismatch in material properties across them is significant, or when reflecting materials are present. Recently Côté et al. reported [8

] a biophotonics code with a triangle based geometry description and polarisation support. This development allows for flexible shape based models, but the lack of spatial hierarchy in the geometry description limits the complexity of the manageable situations.

In this paper we present a Monte Carlo environment suitable for addressing heterogeneous complex tissues in an efficient manner. The geometry is based on a set of triangle meshes, which are structured using a space partitioning scheme in order to improve computational efficiency. Each triangle has information on the two materials immediately in contact with it on either side of its normal, and the combination of these elements into a mesh is used to define a complete interface between the two media. This shape based representation is more flexible than voxel based approaches and allows for more accurate modelling of reflection and refraction caused by a surface. The optimisation techniques applied make it possible to track the paths of large numbers of photons, since the computational time of the photon-geometry intersection test decreases from asymptotically proportional to the number of triangles to logarithmic with the number of triangles.

Previous numerical work on tissue heterogeneity on scales under the average tissue diffusion length has mostly focused on the effect of blood vessels, mainly due to their pervasiveness in body tissues and their strong contrast in absorption relative to surrounding media. Due to the fact that the diffusion approximation is not applicable, Monte Carlo has been the method of choice. The evaluation of the impact of tissue structure on optical diagnostic methods [16

A. Talsma, B. Chance, and R. Graaff, “Corrections for inhomogeneities in biological tissue caused by blood vessels,” J. Opt. Soc. Am. A 18, 932–939 (2001). [CrossRef]

, 17

R. van Veen, W. Verkruysse, and H. Sterenborg, “Diffuse-reflectance spectroscopy from 500 to 1060 nm by correction for inhomogeneously distributed absorbers,” Opt. Lett. 27(4), 246–248 (2002). [CrossRef]

] and treatments [18

J. Barton, T. Pfeffer, and A. Welch, “Optical Monte Carlo modeling of a true port wine stain anatomy,” Opt. Express 2(9), 391–396 (1998). [CrossRef] [PubMed]

] shows that non-homogeneity has often a non-negligible impact on average transport properties and radiation effects. The effect of macroscopic heterogeneity has also been investigated by means of Monte Carlo methods. In this case, MC studies have sometimes concentrated on the effect of transparent layers [19

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

, 7

, 20

Y. Fukui, Y. Ajichi, and E. Okada, “Monte Carlo prediction of near-infrared light propagation in realistic adult and neonatal head models,” Appl. Opt. 42(16), 2881–2887 (2003). [CrossRef] [PubMed]

], as these problems can be modelled with Finite Element Methods only with additional assumptions and computational devices [21

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

]. Further research subjects in this category are broad in scope and include modelling of photoacoustic signals from large vessels [22

I. Patrikeev, Y. Petrov, I. Petrova, D. Prough, and R. Esenaliev, “Monte Carlo modeling of optoacoustic signals from human internal jugular veins,” Appl. Opt. , 46(21), 4820–4827 (2007). [CrossRef] [PubMed]

] and calculations of fluence rate distribution for Photodynamic Therapy around human airways [23

L. Murrer, H. Marijnissen, and W. Start, “Monte Carlo simulations for endobronchial photodynamic therapy: the influence of variations in optical and geometrical properties and of realistic and eccentric light sources,” Laser. Surg. Med. 22, 193–206 (1998). [CrossRef]

We will target trabecular bone as an example of the kind of complex heterogeneous media that can be considered for simulation, demonstrating the code’s possibilities. The versatility of the geometry description presented can be applied to the problems mentioned in the previous paragraph, helping understand the impact of tissue morphological complexity in the field of biomedical optics. However, for evaluation purposes, the characteristic structure of trabecular bone offers an excellent test bench for the program’s capabilities.

Despite the illustrative purposes of the simulation results presented, the ability to study light propagation in bone also has importance in itself. Optical investigation of osseous tissue has received increased attention in the last years [24

A. Pifferi, A. Torricelli, P. Taroni, A. Bassi, E. Chikoidze, E. Giambattistelli, and R. Cubeddu, “Optical biopsy of bone tissue: a step toward the diagnosis of bone pathologies,” J. Biomed. Opt. 9(3), 474–479 (2004). [CrossRef] [PubMed]

]. This interest is probably due to the inherent low invasiveness, potential low cost and intrinsic high information content, derived from the characteristic spectra of individual tissue components. Some diseased states of bone, like osteoporosis, are excellent candidates for optical diagnosis. Osteoporosis is defined [25

M. Beers and R. Berkow, eds. The Merck Manual of Diagnosis and Therapy , 17th ed. (Merck Research Laboratories, Whitehouse Station, NJ, 1999).

] as a generalised, progressive diminution of bone density (bone mass per unit volume), causing skeletal weakness. Its detection is currently achievable by means of biopsy, serum analysis and radiographic techniques. Due to some drawbacks of the mentioned methods, including invasive-ness, high cost, low portability and high radiation dose, interest in novel methods for diagnosis of diseased states of bone is on the rise. In recent times, quantitative ultrasound has emerged as a potential diagnosis method [26

F. Hartl, A. Tyndall, M. Kraenzlin, C. Bachmeier, C. Gückel, U. Senn, D. Hans, and R. Theiler, “Discriminatory Ability of Quantitative Ultrasound Parameters and Bone Mineral Density in a Population-Based Sample of Post-menopausal Women With Vertebral Fractures: Results of the Basel Osteoporosis Study,” J. Bone Miner. Res. 17, 321–330 (2002). [CrossRef] [PubMed]

], though there is still some debate in the literature about its actual value for screening purposes [27

S. Nayak, I. Olkin, H. Liu, M. Grabe, M. Gould, I. Allen, D. Owens, and D. Bravata, “Meta-Analysis: Accuracy of Quantitative Ultrasound for Identifying Patients with Osteoporosis,” Ann. Intern. Med. 144(11), 832–841 (2006). [PubMed]

In this line, some in vitro [28

A. Takeuchi, R. Araki, S. Proskurin, Y. Takahashi, Y. Yamada, J. Ishii, S. Katayama, and A. Itabashi, “A New Method of Bone Tissue Measurement Based upon Light Scattering,” J. Bone Miner. Res. 12(2), 261–266 (1997). [CrossRef] [PubMed]

] and in vivo [24

] studies have tried to relate the scattering and absorption properties of bone in the Near Infrared (NIR) to its quality. Though these studies use very promising methods, conclusive results have not yet been obtained. Despite the sophisticated experimental methods applied, these reports rely on quite simple mechanisms of light transport within bone. They generally assume that this tissue can be regarded as a homogeneous isotropic turbid medium and that its mean optical properties (especially absorption) can be constructed as a weighted average from those of its basic constituents. It is also implicitly assumed that the weighing of the individual spectra is a direct measure of their relative concentration in tissue.

This is in fact a very common premise in macroscopic models of tissue, and it implicitly implies a homogeneous distribution of scattering and absorption centres at dimensions below the diffusion length. As mentioned before, much of previous numerical work on tissue heterogeneity has focused on this problem for the case of blood vessels. In fact, the complex morphology of trabecular bone, whose complexity clearly departs from the assumed homogeneity principle also calls for further investigation. The effect of this microscopic structure on the behaviour of light under the diffusion limit is unclear, though some studies on other heterogeneous media with structural regularity [30

J. Hebden, J.J. García Guerrero, V. Chernomordik, and A.H. Gandjbakhche, “Experimental evaluation of an anisotropic scattering model of a slab geometry,” Opt. Lett. 29, 2518–2520 (2004). [CrossRef] [PubMed]

], suggest that it could be significant. Here, we show how the presented code can be used to efficiently compute light transport in cancellous bone based on models reconstructed from micro Computer Tomography (μCT) data. This reality-based computation of photon migration in a complex two-phase heterogeneous media would not have been possible to the same level of detail without the advances presented in this paper.

In the following, the principles on which the implementation is based are first introduced. The steps taken for validating the environment are discussed next, together with some performance measures that show the impact of problem geometry on computational time. As an example of the possibilities of the new environment, simulations of light transport in a model of trabecular bone as recovered from μCT scans are presented.

2. Implementation

The fundamentals of Monte Carlo simulation of light transport in turbid media have been discussed extensively in the literature. This implementation is based on the work by Wang et al. [5

L. Wang, S. Jacques, and L. Zheng, “MCML - Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Meth. Prog. Bio. 47, 131–146 (1995). [CrossRef]

] for multilayered geometries, with an extension to arbitrary 3D structures based on collections of triangles. In this respect, our environment is similar to the one described by Côté and Vitkin. However, the organisation of the geometry elements in a bounding volumes hierarchy allows computation of very large meshes. The application of this strategy is to the best of our knowledge novel in the field of biomedical optics.

The structure of the environment architecture has been depicted in Fig. 1. It is based on a MC core coded in the C programming language, which interfaces with the pre- and post-processing steps through files based on the Hierarchical Data Format HDF5. Usage of this standard data organisation scheme [31

“HDF5 - A New Generation of HDF,” Available at http://hdf.ncsa.uiuc.edu/HDF5 (2000). NCSA, University of Illinois at Urbana Champaign.

] allows an efficient implementation in terms of I/O and ensures compatibility with widespread tools and libraries. This encapsulation of the MC core gives increased flexibility in the pre- and post-processing steps, which can be in this way performed using more convenient tools, and allowed us to concentrate on an efficient photon path computation. The complete source code of the implementation of the algorithm described in this article is available on the World Wide Web under the name TriMC3D [32

E. Margallo-Balbás, “Home page of TriMC3D,” (2007), http://trimc3d.et.tudelft.nl/

The environment includes a parallel implementation of the photon tracing algorithm based on the Messaging Passing Interface standard, which permits a transparent distribution of a single problem within a computer cluster. This is important for the Monte Carlo technique, which relies on high number of photon paths being calculated for a good signal-to-noise ratio. Although the HDF5 standard includes mechanisms for parallel data access, a collector process taking care of file output was chosen. This option does not require an underlying parallel file system to ensure data consistency and it is appropriate for the moderate data throughput produced by the MC simulations. It was found that the overhead of the parallel implementation was negligible for moderately large numbers of photons.

Fig. 1. Environment architecture for Monte Carlo simulation

Materials are described by means of their scattering and absorption coefficients, together with the anisotropy coefficient of the phase function. Only the common Henyey-Greenstein phase function has been implemented, but the extension to other distributions is straightforward. The geometry description is shape based, consisting in collection of triangles defining the interfaces between different media. In this sense, care has to be taken so as to ensure topological consistency in the model.

Triangles are organised in a hierarchy of bounding volumes, implementing the space partitioning scheme already mentioned. Figure 2 illustrates this concept in the two-dimensional case. Instead of performing individual checks against all elements that make up the geometry, they are organised into a tree structure by repeatedly subdividing the occupied space. In this tree, each node strictly contains all of its siblings, making it unnecessary to further descend the hierarchy if the upper node is known not to intersect the photon step.

The bounding volume hierarchy is based on axis-aligned bounding boxes and it is generated in a pre-processing step. It consists of a tree data structure, in which each node potentially contains both triangles and links to sibbling nodes. During loading of the geometry, a bounding box is pre-computed for each of the tree nodes in a bottom-up manner making use of a recursive function. All nodes are formally equal, with the difference that leaf (terminal) nodes do not have pointers to any further sibblings, only containing a list of associated triangles.

Since general hierarchies with a multi-level distribution of triangles in groups can be processed by the MC core, various axis-aligned space partitioning schemes can be in principle used. Given the volumetric uniformity of the triangle distribution in the problems discussed later in this paper, a “loose” octree organisation [34

A. Glassner, “Space Subdivision for Fast Ray Tracing,” IEEE Comput. Graph. 4(10), 15–22 (1984).

] was adopted. Octrees recursively divide each cubic node into eight equal sub-volumes, being a three-dimensional analogue of the quadtree bidimensional spatial decomposition used for illustration purposes in Fig. 2. A more general space partitioning scheme is the three-dimensional kd-tree, in which space is recursively subdivided by planes normal to the Cartesian axes. The performance of the different space partitioning schemes is application specific and the reader is referred to the literature [15

] for further information.

Fig. 2. Illustration of the concept of space partitioning. The final geometry elements on the left of the geometry are organised into a three-layer (A,B,C) hierarchy that can be efficiently traversed, limiting the number of final photon-element tests to a minimum.

On initialisation of the geometry intersection routine, a bounding box is computed for the maximum possible photon segment, corresponding to undisturbed propagation in the current medium to the full photon step length. The traversal of the tree is carried out from the root node to the leaf elements in a recursive way. At each node, all associated triangles are first checked for intersection with the photon path. In our implementation, this intersection test is composed of a simple box-box overlap check, followed by a ray-triangle intersection evaluation. The photon segment bounding box and the triangle bounding box are evaluated in the first stage. The ray-triangle intersection test has been adapted from the algorithm proposed by Möller et al. [33

T. Möller and B. Trumbore, “Fast, minimum storage ray-triangle intersection,” J. Graph. Tools 2(1), 21–28 (1997).

], which offers a fast result while not requiring storage of pre-computed coefficients. If a triangle is found to intersect the photon path, it is recorded as the closest point in the geometry in the current photon direction and its distance to the current photon position is used to recompute the maximum possible photon segment and its bounding box. After this, the routine resumes checking the remaining triangles in the node, if any.

After triangles, sibbling nodes in the hierarchy of bounding volumes are evaluated for intersection. This is done through a fast box-box overlap test, followed by a more restrictive ray-box intersection check, adapted from the implementation proposed by Williams [35

A. William, S. Barrus, R. Morley, and P. Shirley, “An Efficient and Robust Ray-Box Intersection Algorithm,” J. Graph. Tools 10(1), 49–54 (2005). [CrossRef]

]. It is interesting that for many applications in biomedical optics, the photon mean free path will be much smaller that the minimum element size. Under these circumstances, a direct ray-box intersection test results in a very inefficient tree traversal. This contrasts with the more common situation in computer graphics, where ray length is often assumed to be comparable to the geometry size or even infinite. We have found that the described test sequence maintains a uniform performance independently of the size of the photon step relative to the element size.

If no intersection is found, propagation of the photon is performed following the MCML algorithm [5

L. Wang, S. Jacques, and L. Zheng, “MCML - Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Meth. Prog. Bio. 47, 131–146 (1995). [CrossRef]

]. If the result of the check is positive, the closest triangle and the minimum distance to intersection (d) are returned. The photon is then propagated to the interface and the remaining photon step length is decreased by dμ_t , where μ_t = μ_a +μ_s is the interaction coefficient of the current medium. After this, a choice is made as whether reflection or transmission takes place. This follows Wang’s [5

L. Wang, S. Jacques, and L. Zheng, “MCML - Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Meth. Prog. Bio. 47, 131–146 (1995). [CrossRef]

] approach, in which the reflection coefficient is compared to a uniform random number and the full weight of the photon is either reflected or transmitted. The reflection coefficient (R) for an incident angle α_i relative to the triangle normal is 1 for α_i > α_cr , where α_cr is the critical angle. For smaller angles, R is obtained from Fresnel’s formulas, as shown in Eq. 1, where polarisation is not taken into account. In Eq. 1, α_t is the angle of transmission, as given by Snell’s law. The photon direction is updated depending on whether it should be reflected or refracted according to the relations of ray optics.

Fig. 3. Separate validation of the ray-tracing algorithm and the photon transport code. (a) Spherical convergent lens. (b) Normalised reflectance as a function of the polar angle for a homogeneous slab, compared to MCML results.

R (α_{i}) = \frac{1}{2} [\frac{\sin^{2} (α_{i} - α_{t})}{\sin^{2} (α_{i} + α_{t})} + \frac{\tan^{2} (α_{i} - α_{t})}{\tan^{2} (α_{i} + α_{t})}]

(1)

Once a photon-surface intersection is processed, the affected triangle is excluded from the geometry for the following check. This is done in order to prevent faulty results due to floating point errors, as after updating its position the photon lies directly on the triangle plane. Once a scattering event takes place, or another reflection/refraction happens, the restriction on this triangle is removed.

Since the environment is meant for arbitrary geometries, no simplifications can be made based on the problem spatial invariance and the convolutional schemes for extended sources described by Wang et al.[36

L. Wang, S. Jacques, and L. Zheng, “CONV - convolution for responses to a finite diameter photon beam incident on multi-layered tissues,” Comput. Meth. Prog. Bio. 54, 141–150 (1997). [CrossRef]

] cannot be used. To address this issue, the environment is able to directly model Gaussian and isotropic sources, besides pencil beam excitation. Extension to other emission patterns is straightforward. The modelling of Gaussian beams has been adopted from Tycho’s work [37

A. Tycho, T. Jørgensen, H. Yura, and P. Andersen, “Derivation of a Monte Carlo method for modeling heterodyne detection in optical Coherence tomography systems,” Appl. Opt. 41(31), 6676–6691 (2002). [CrossRef] [PubMed]

]. Ray launching times are in most situations negligible compared to the cost of photon path calculation, even for the more sophisticated source radiation patterns. However, launching angles and positions have to be explicitly sampled and the filtering effect of the convolution is lost, which often translates into a need for larger number of simulated photon paths. This is in fact caused by the lack of space invariance in the geometry more than due to limitations in the code. Detectors are modelled explicitly in the geometry as exit surfaces on which photon paths are recorded. Several parameters of interest are registered, such as the photon exit coordinates, the cosine of its angle with respect to the surface’s normal, the path-length from the moment of launching, the accumulated time of flight and the remaining weight.

3. Validation and performance

The basic validity of light transport computation based on Monte Carlo methods and the physics of light propagation across medium interfaces are well-established. Therefore, the 3D ray tracing algorithm has been validated separately from the Monte Carlo routines of absorption and scattering. Numerous well-known problems of geometrical optics have been used to check the reliability of the surface intersection test and the accuracy of the calculated intensity and direction of the refracted and reflected rays. The convergence of a beam of parallel rays incident on a lens composed of 2.9×10⁴ elements to its paraxial focal point is depicted in Fig. 3(a), where the problem geometry has been overlaid to the photon paths. It can be observed that the effect of spherical aberration is reproduced for rays far from the optical axis. The source intensities have been individually validated against their analytical models. Good agreement was found for sufficiently high numbers of photons.

Fig. 4. Normalised computational time as a function of model complexity, showing the effect of hierarchical space partitioning.

The algorithms modelling light transport in scattering media have been validated against the results of the commonly used MCML code [5

L. Wang, S. Jacques, and L. Zheng, “MCML - Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Meth. Prog. Bio. 47, 131–146 (1995). [CrossRef]

]. These results have been compared to the diffusion equation and have been used as a reference for other validations in the past [8

]. Figure 3(b) shows the angle resolved reflectance generated by a infinitely narrow source incident on a homogeneous slab of 0.2mm thickness with μ_s = 9mm ^-1, μ_a = 1mm ^-1 and g = 0.75. Both solutions were computed with 5×10⁶ photons and show an excellent agreement, with the residual error due to the still finite number of simulated paths. This conformity is not surprising, as the basic algorithms for scattering and absorption were implemented on the basis of Wang’s [5

L. Wang, S. Jacques, and L. Zheng, “MCML - Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Meth. Prog. Bio. 47, 131–146 (1995). [CrossRef]

] article.

In order to evaluate the performance of the presented Monte Carlo code, a simple model of a spherical absorber in an otherwise homogeneous slab was taken as reference problem. The sphere has the optical properties μ_s = 10mm ^-1, μ_a = 0.01mm ^-1 and an anisotropy coefficient of g = 0.9, while the surrounding medium can be characterised by the coefficients μ_s = 10mm ^-1, μ_a = 0.05mm ^-1 and an anisotropy factor g = 0.9. Figure 4 represents the computational time as a function of the number of triangles used to approximate the sphere. For ease of comparison, the simulation times were normalised to the time needed for a 76 triangle model.

In the considered scenario, and making use of an octree hierarchical space organisation, the photon-geometry intersection test requires negligible time when compared to the rest of the photon tracking functions. This effect is reflected in the fact that no observable increase in the computation time is produced by the growing complexity. This can be achieved because only few costly ray-triangle intersections need to be computed after the tree traversal for each photon step, independent of the actual level of detail in the sphere description. We have found that this effect can be found in most practical situations, with the geometry being in practice neutral with respect to computational time.

For comparison purposes, Fig. 4 also shows the evolution of computational time for the spherical absorbing inclusion model if a flat (non-hierarchical) organisation of triangles is applied. It becomes apparent that already for moderate levels of model complexity (~ 10³) the computational cost introduced by the photon-geometry intersection test becomes prohibitive. In order to observe the logarithmic dependency in computational time due to the hierarchical organisation alone, the time needed to conduct a photon-surface intersection test has been isolated. The circles in Fig. 4 represent the normalised time as a function of the model complexity. The curve clearly reproduces the expected behaviour, approximately following the logarithmic dependency represented by the fit to a straight line on the lin-log graph.

4. Application to trabecular bone

Triangle meshes reconstructed from μCT data have been simulated. Bone samples were obtained from the trabecular space of vertebrae and from the femoral head of adult Dutch milk goats. The animals were sacrificed after a long term study. Remaining tissues had been frozen after death and specimens were retrieved and cut to a suitable size for analysis. Three-dimensional data was obtained using a μCT system (μCT20, Scanco Medical AG, Zürich, Switzerland). Briefly, the samples were placed in a fluid to avoid dehydration and scanned at a resolution of 18μm [38

E. Giesen and T. van Eijden, “The three-dimensional cancellous bone architecture of the human mandibular condyle,” J. Dent. Res. , 79, 957–963 (2000). [CrossRef] [PubMed]

] Spherical volumes of interest were extracted from the voxel-based dataset, and were then thresholded and transformed into a triangular mesh using a marching cubes algorithm. For a sphere with a diameter of 7.5mm, approximately 7 million triangles were generated after this step. In all cases, the resulting mesh was then simplified to 10% of its original number of triangles for later processing convenience using the quadric decimation metric first described by Garland et al. [39

M. Garland and P. Heckbert, “Surface Simplification Using Quadric Error Metrics,” in SIGGRAPH (1997).

], without noticeable degradation of the resulting model.

For this report we assume that the calcified matrix has similar optical properties to those of cortical bone, with n = 1.6, μ_a = 0.1mm ^-1, g = 0.93 and μ_s = 27mm ^-1. The assumed values are representative for an excitation wavelength around 800nm according to Firbank [40

M. Firbank, M. Hiraoka, M. Essenpreis, and D. Delpy, “Measurement of the optical properties of the skull in the wavelength range 650–950 nm,” Phys. Med. Biol. 38, 503–510 (1993). [CrossRef] [PubMed]

], among others [29

N. Ugryumova, S. Matcher, and D. Attenburrow, “Measurement of bone mineral density via light scattering,” Phys. Med. Biol. 49, 469–483 (2004). [CrossRef] [PubMed]

, 24

Table 1. Assumed optical properties for the two constituents of bone

Material	n	μ_s	μ_a	g
Calcified matrix	1.6	27mm ^-1	0.025mm ^-1	0.93
Marrow (40°C)	1.4	1.5mm ^-1	0.1mm ^-1	0.95

The inner space is filled with marrow, whose properties have not been experimentally determined to date. For the purposes of illustration, they will be assumed similar to those of adipose tissue in other parts of the body. Bashkatov [41

G. Bashkatov, “Optical properties of human skin, subcutaenous and mucous tissues in the wavelength range from 400 to 2000nm,” J. Phys. D Appl. Phys. 38, 2543–2555 (2005). [CrossRef]

] has reported a relatively high scattering coefficient (μ′ _s = 1mm ^-1) at λ = 800nm for human subcutaneous adipose tissue at room temperature (T ~ 20°C). In contrast with this value, Tsai [42

C. Tsai, Y. Yang, C. Han, J. Hsieh, and M. Chang, “Measurement and simulation of light distribution in biological tissues,” Appl. Opt. 40(31), 5770–5777 (2001). [CrossRef]

] and van Veen [43

R. van Veen, H. Sterenborg, A. Pifferi, A. Torricelli, E. Chikoidze, and R. Cubeddu, “Determination of visible near-IR absorption coefficients of mammalian fat using time and spatially resolved diffuse reflectance and transmission spectroscopy,” J. Biomed. Opt. 10(5) (2005). [PubMed]

] have shown that for lipids alone this coefficient may decrease significantly when temperatures are closer to body conditions. Tsai’s estimations of the scattering coefficient and the anisotropy factor of porcine fat at 40°C are μ_s = 1.5mm ^-1,g = 0.95. Though temperature may play a role in the differences in reported results, it is important to note that physiological adipose tissue contains several other components contributing to scattering, including vessels, nerves and collagen fibres. Further, lipids in the body are in a different state from ex-vivo processed samples. How the actual state of lipids in interstitial marrow may condition its scattering coefficient is uncertain. Absorption properties seem to have a lower temperature dependency [43

] and for the rest we will use the value reported by Bashkatov (μ_a = 0.1mm ^-1) for this coefficient. The issue of the optical properties of marrow is further complicated because other non-considered tissue constituents may play an important role and two types of marrow (red and yellow) with strongly varying composition exist. For the purposes of this article, i.e. illustration of the possibilities of the presented code, the question of exact optical properties plays a secondary role. Table 1 summarises the set of assumed transport coefficients that were used in this study.

Fig. 5. Simulation of light transport in trabecular bone

Figure 5(a) shows the model of a reconstructed bone sphere with d = 5mm together with the corresponding simulation geometry. The interstitial medium (marrow) is in this case supposed to fill the space. The fluence distribution was obtained from 8 × 10⁶ photons launched from a pencil-beam source in the lower part of the figure. The values represented in the figures are the decimal logarithm of the number of photons scored at each grid element. Contour lines with constant fluence have been added to the figures, while the thinner dotted lines represent the limits of the trabeculae. Though the volume of interest is not very big, the influence of trabeculae on the long-range distribution of fluence in tissue can be seen. The large number of photons calculated provides a good signal-to-noise ratio and high-quality fluence maps. These simulations were computed within 1 hour of CPU time on a 3.2GHz Pentium 4 platform.

Time resolved curves through a bone sphere with a diameter d = 5mm have been computed. A Gaussian source with a diameter at 1/e ² intensity of 0.5mm in the launching plane was used in this case for excitation, with its focus set at infinite. The detector was placed on the opposite side of the sphere relative to its centre. A 0.5mm margin was left between the sphere and the collecting and emitting elements. Detectors had a diameter of 0.7mm and were assumed to collect light incident from all directions. This assumed lack of numerical aperture could be reasonable for high acceptance angle detectors, such as photodiodes, but it could be unrealistic for fibre based systems. Within these simulations, we found that for NA=0.48 and in a medium with refractive index n = 1.4, the number of simulated photons should be increased by a factor of approximately 5 to get the same accumulated photon weight on the detector. If very small numerical apertures and detectors have to be taken into account, a different approach using the reciprocity theorem of the RTE should be considered, as direct application of Monte Carlo wastes a large number of photons. This is not specific of the presented code, but rather characteristic of one-pass Monte Carlo methods.

The source-detector pair was rotated in a plane around the sphere in 10° increments and time resolved curves with ~ 4×10⁷ photons were acquired for each position. Due to the large number of paths needed for good temporal resolution and the numerous sources involved in the simulation, the computation was performed in the nodes of the DAS2 parallel supercomputer [44

H. Bal, R. Bhoedjang, R. Hofman, C. Jacobs, T. Kielmann, J. Maassen, R. van Nieuwpoort, J. Romein, L. Renambot, T. Rühl, R. Veldema, K. Verstoep, A. Baggio, G. Ballintijn, I. Kuz, G. Pierre, M. van Steen, A. Tanenbaum, G. Doornbos, D. Germans, H. Spoelder, E.-J. Baerends, S. van Gisbergen, H. Afsermanesh, D. van Albada, A. Belloum, D. Dubbeldam, Z. Hendrikse, B. Hertzberger, A. Hoekstra, K. Iskra, D. Kandhai, D. Koelma, F. van der Linden, B. Overeinder, P. Sloot, P. Spinnato, D. Epema, A. van Gemund, P. Jonker, A. Radulescu, C. van Reeuwijk, H. Sips, P. Knijnenburg, M. Lew, F. Sluiter, L. Wolters, H. Blom, C. de Laat, and A. van der Steen, “The distributed ASCI Supercomputer project,” SIGOPS Oper. Syst. Rev. 34(4), 76–96 (2000). [CrossRef]

], with a total number of 1.5 × 10⁹ photons traced in ~ 12h. Figure 6 shows some curves representative of the obtained set of transient responses. It is unclear at this point how the angular dependency found in the curves relates to the underlying model properties. More accurate investigation on this issue is needed, but the example application presented successfully illustrates the potential of the presented code.

Fig. 6. Simulated impulse responses across a sphere of trabecular bone for different relative angles

5. Conclusions

In this paper we present a new Monte Carlo code based on a shape based geometry description. Drawing from concepts from computer graphics and previous work on stochastic modelling of radiative transport in tissue optics, we have been able to address complex geometries in a general, flexible and efficient way. We have found that a hierarchical organisation of the triangle mesh, together with a suitable intersection test strategy delivers a good computational efficiency. This kind of cross-seeding between different disciplines making use of Monte Carlo methods for computation of radiative transfer has been already highlighted in the literature [13

E. Veach, “Robust Monte Carlo methods for light transport simulation,” Ph.D. thesis, Stanford University (1997).

]. A complete validation of the code has been conducted and the main results have been presented.

Upper bounds to the complexity of the geometry in space partitioning schemes like the one proposed come from memory limitations more than from exploding computational cost. Effort has been devoted to efficient memory usage in the implementation, but an upper limit exists, as the whole model should fit in the main computer memory to avoid page swapping to disk. With current computers these limits are far even for complex models like the ones reconstructed from imaging devices. Using a 32 bit memory space, at least 60 million triangles are tractable. If larger meshes are required, memory management techniques taking spatial relevance into account could be introduced. A departure from the standard sequential photon computation scheme [45

M. Pharr, C. Kolb, R. Gershbein, and P. Hanrahan, “Rendering Complex Scenes with Memory-Coherent Ray Tracing,” in SIGGRAPH ’97: Proceedings of the 24th annual conference on Computer graphics and interactive techniques (ACM Press/Addison-Wesley Publishing Co., 1997), pp. 101–108. [CrossRef]

] can also lead to an improvement in memory locality.

In this implementation the geometry and material properties are assumed to remain constant during light propagation in tissue. This means that an assumption of temporal invariance underlies the current implementation. Sources with a finite temporal response function can be easily considered by means of time-domain convolution. However, ultrafast events occurring in the geometry or material properties due to the light pulse propagation cannot be easily modelled using the described environment. Although this limitation is not important for standard tissue spectroscopy techniques and many other applications, an extension to include time-dependent geometry could be useful for ultrafast studies.

As a relevant and demanding application, we have investigated the effect of the micro structure of trabecular bone on its long-range transport properties. Direct simulation on bone models reconstructed from μCT data was successfully conducted and spatial and temporal information about light propagation in tissue was obtained. Lack of sound experimental information regarding the properties of the individual components of bone currently limits a more detailed numerical investigation. Future work will also include direct comparison of simulation results to time-resolved measurements on the modelled sample.

Acknowledgements

We gratefully acknowledge Dr. Th. Smit and Mr. A.J. Veen from the Department of Physics and Medical Technology at the Vrije Universiteit Medical Centre in Amsterdam for kindly providing us with bone specimens. We are also thankful to Mr. L.J. van Ruijven from the Department of Functional Anatomy at the Amsterdam Medical Centre for his work on the acquisition of the μCT scans used in this study. We want to express our gratitude to Dr. Antonio Pifferi and Dr. Paola Taroni for encouraging and instructive discussions. This work has been supported by the Dutch Technology Foundation, STW, Applied Science Foundation of NOW and the technology programme of the Dutch Ministry of Economic Affairs under contract number STW7505 and by the Delft Centre for Mechatronics and Microsystems.

References and links

1.	X. Gu, Q. Zhang, M. Bartlett, L. Schutz, L. Fajardo, and H. Jiang, “Differentiation of cysts from solid tumours in the breast with Diffuse Optical Tomography,” Acad. Radiol. 11, 53–60 (2004). [CrossRef] [PubMed]
2.	P. Taroni, G. Danesini, A. Torricelli, A. Pifferi, L. Spinelli, and R. Cubbedu, “Clinical trial of time-resolved scanning optical mammography at 4 wavelengths between 683 and 975 nm,” J. Biomed. Opt. 9(3), 464–473 (2004). [CrossRef] [PubMed]
3.	S. Brown, E. Brown, and I. Walker, “The present and future role of photodynamic therapy in cancer treatment,” Oncology 5, 497–508 (2004). [PubMed]
4.	S. Chandrasekhar, Radiative Transfer (Oxford University Press, London, 1960).
5.	L. Wang, S. Jacques, and L. Zheng, “MCML - Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Meth. Prog. Bio. 47, 131–146 (1995). [CrossRef]
6.	T. Pfefer, J. Barton, E. Chan, M. Ducros, B. Sorg, T. Milner, J. Nelson, and A. Welch, “A three-dimensional modular adaptable grid numerical model for light propagation during laser irradiation of skin tissue,” IEEE J. Sel. Top. Quantum Electron. 2(4), 934–942 (1996). [CrossRef]
7.	D. Boas, J. Culver, J. Stott, and A. Dunn, “Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head,” Opt. Express 10(3), 159–170 (2002). [PubMed]
8.	D. Côté and I. Vitkin, “Robust concentration determination of optically active molecules in turbid media with validated three-dimensional polarization sensitive Monte Carlo calculations,” Opt. Express 13, 148–163 (2005). [CrossRef] [PubMed]
9.	F. Brown, R. Barrett, T. Booth, J. Bull, L. Cox, R. Forster, J. Goorley, R. Mosteller, S. Post, R. Prael, E. Selcow, A. Sood, and J. Sweezy, “MCNP Version 5,” Trans. Am. Nucl. Soc. 87, 273 (2002).
10.	J. Howell, “The Monte Carlo method in radiative heat transfer,” J. Heat Transf. 120(3), 547–560 (1998). [CrossRef]
11.	A. Marshak and A. B. Davis, eds., 3D Radiative Transfer in Cloudy Atmospheres (Springer-Verlag, Berlin, 2005). [CrossRef]
12.	A. S. Glassner, Principles of Digital Image Synthesis (Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 1994).
13.	E. Veach, “Robust Monte Carlo methods for light transport simulation,” Ph.D. thesis, Stanford University (1997).
14.	E. Cerezo, F. Pérez, X. Pueyo, F. Serón, and F. Sillion, “A survey on participating media rendering techniques,” Visual. Comput. 21, 303–328 (2005). [CrossRef]
15.	L. Szirmay-Kalos, V. Havran, B. Balász, and L. Szécsi, “On the efficiency of ray-shooting acceleration schemes,” in SCCG ’02: Proceedings of the 18th spring conference on Computer graphics , pp. 97–106 (ACM Press, New York, NY, USA, 2002). [CrossRef]
16.	A. Talsma, B. Chance, and R. Graaff, “Corrections for inhomogeneities in biological tissue caused by blood vessels,” J. Opt. Soc. Am. A 18, 932–939 (2001). [CrossRef]
17.	R. van Veen, W. Verkruysse, and H. Sterenborg, “Diffuse-reflectance spectroscopy from 500 to 1060 nm by correction for inhomogeneously distributed absorbers,” Opt. Lett. 27(4), 246–248 (2002). [CrossRef]
18.	J. Barton, T. Pfeffer, and A. Welch, “Optical Monte Carlo modeling of a true port wine stain anatomy,” Opt. Express 2(9), 391–396 (1998). [CrossRef] [PubMed]
19.	E. Okada, M. Firbank, M. Schweiger, S. Arridge, M. Cope, and D. Delpy, “Theoretical and experimental investigation of near-infrared light propagation in a model of the adult head,” Appl. Opt. 36(1), 21–31 (1997). [CrossRef] [PubMed]
20.	Y. Fukui, Y. Ajichi, and E. Okada, “Monte Carlo prediction of near-infrared light propagation in realistic adult and neonatal head models,” Appl. Opt. 42(16), 2881–2887 (2003). [CrossRef] [PubMed]
21.	M. Firbank, S. Arridge, M. Schweiger, and D. Delpy, “An investigation of light transport through scattering bodies with non-scattering regions,” Phys. Med. Biol. 41, 767–783 (1996). [CrossRef] [PubMed]
22.	I. Patrikeev, Y. Petrov, I. Petrova, D. Prough, and R. Esenaliev, “Monte Carlo modeling of optoacoustic signals from human internal jugular veins,” Appl. Opt. , 46(21), 4820–4827 (2007). [CrossRef] [PubMed]
23.	L. Murrer, H. Marijnissen, and W. Start, “Monte Carlo simulations for endobronchial photodynamic therapy: the influence of variations in optical and geometrical properties and of realistic and eccentric light sources,” Laser. Surg. Med. 22, 193–206 (1998). [CrossRef]
24.	A. Pifferi, A. Torricelli, P. Taroni, A. Bassi, E. Chikoidze, E. Giambattistelli, and R. Cubeddu, “Optical biopsy of bone tissue: a step toward the diagnosis of bone pathologies,” J. Biomed. Opt. 9(3), 474–479 (2004). [CrossRef] [PubMed]
25.	M. Beers and R. Berkow, eds. The Merck Manual of Diagnosis and Therapy , 17th ed. (Merck Research Laboratories, Whitehouse Station, NJ, 1999).
26.	F. Hartl, A. Tyndall, M. Kraenzlin, C. Bachmeier, C. Gückel, U. Senn, D. Hans, and R. Theiler, “Discriminatory Ability of Quantitative Ultrasound Parameters and Bone Mineral Density in a Population-Based Sample of Post-menopausal Women With Vertebral Fractures: Results of the Basel Osteoporosis Study,” J. Bone Miner. Res. 17, 321–330 (2002). [CrossRef] [PubMed]
27.	S. Nayak, I. Olkin, H. Liu, M. Grabe, M. Gould, I. Allen, D. Owens, and D. Bravata, “Meta-Analysis: Accuracy of Quantitative Ultrasound for Identifying Patients with Osteoporosis,” Ann. Intern. Med. 144(11), 832–841 (2006). [PubMed]
28.	A. Takeuchi, R. Araki, S. Proskurin, Y. Takahashi, Y. Yamada, J. Ishii, S. Katayama, and A. Itabashi, “A New Method of Bone Tissue Measurement Based upon Light Scattering,” J. Bone Miner. Res. 12(2), 261–266 (1997). [CrossRef] [PubMed]
29.	N. Ugryumova, S. Matcher, and D. Attenburrow, “Measurement of bone mineral density via light scattering,” Phys. Med. Biol. 49, 469–483 (2004). [CrossRef] [PubMed]
30.	J. Hebden, J.J. García Guerrero, V. Chernomordik, and A.H. Gandjbakhche, “Experimental evaluation of an anisotropic scattering model of a slab geometry,” Opt. Lett. 29, 2518–2520 (2004). [CrossRef] [PubMed]
31.	“HDF5 - A New Generation of HDF,” Available at http://hdf.ncsa.uiuc.edu/HDF5 (2000). NCSA, University of Illinois at Urbana Champaign.
32.	E. Margallo-Balbás, “Home page of TriMC3D,” (2007), http://trimc3d.et.tudelft.nl/
33.	T. Möller and B. Trumbore, “Fast, minimum storage ray-triangle intersection,” J. Graph. Tools 2(1), 21–28 (1997).
34.	A. Glassner, “Space Subdivision for Fast Ray Tracing,” IEEE Comput. Graph. 4(10), 15–22 (1984).
35.	A. William, S. Barrus, R. Morley, and P. Shirley, “An Efficient and Robust Ray-Box Intersection Algorithm,” J. Graph. Tools 10(1), 49–54 (2005). [CrossRef]
36.	L. Wang, S. Jacques, and L. Zheng, “CONV - convolution for responses to a finite diameter photon beam incident on multi-layered tissues,” Comput. Meth. Prog. Bio. 54, 141–150 (1997). [CrossRef]
37.	A. Tycho, T. Jørgensen, H. Yura, and P. Andersen, “Derivation of a Monte Carlo method for modeling heterodyne detection in optical Coherence tomography systems,” Appl. Opt. 41(31), 6676–6691 (2002). [CrossRef] [PubMed]
38.	E. Giesen and T. van Eijden, “The three-dimensional cancellous bone architecture of the human mandibular condyle,” J. Dent. Res. , 79, 957–963 (2000). [CrossRef] [PubMed]
39.	M. Garland and P. Heckbert, “Surface Simplification Using Quadric Error Metrics,” in SIGGRAPH (1997).
40.	M. Firbank, M. Hiraoka, M. Essenpreis, and D. Delpy, “Measurement of the optical properties of the skull in the wavelength range 650–950 nm,” Phys. Med. Biol. 38, 503–510 (1993). [CrossRef] [PubMed]
41.	G. Bashkatov, “Optical properties of human skin, subcutaenous and mucous tissues in the wavelength range from 400 to 2000nm,” J. Phys. D Appl. Phys. 38, 2543–2555 (2005). [CrossRef]
42.	C. Tsai, Y. Yang, C. Han, J. Hsieh, and M. Chang, “Measurement and simulation of light distribution in biological tissues,” Appl. Opt. 40(31), 5770–5777 (2001). [CrossRef]
43.	R. van Veen, H. Sterenborg, A. Pifferi, A. Torricelli, E. Chikoidze, and R. Cubeddu, “Determination of visible near-IR absorption coefficients of mammalian fat using time and spatially resolved diffuse reflectance and transmission spectroscopy,” J. Biomed. Opt. 10(5) (2005). [PubMed]
44.	H. Bal, R. Bhoedjang, R. Hofman, C. Jacobs, T. Kielmann, J. Maassen, R. van Nieuwpoort, J. Romein, L. Renambot, T. Rühl, R. Veldema, K. Verstoep, A. Baggio, G. Ballintijn, I. Kuz, G. Pierre, M. van Steen, A. Tanenbaum, G. Doornbos, D. Germans, H. Spoelder, E.-J. Baerends, S. van Gisbergen, H. Afsermanesh, D. van Albada, A. Belloum, D. Dubbeldam, Z. Hendrikse, B. Hertzberger, A. Hoekstra, K. Iskra, D. Kandhai, D. Koelma, F. van der Linden, B. Overeinder, P. Sloot, P. Spinnato, D. Epema, A. van Gemund, P. Jonker, A. Radulescu, C. van Reeuwijk, H. Sips, P. Knijnenburg, M. Lew, F. Sluiter, L. Wolters, H. Blom, C. de Laat, and A. van der Steen, “The distributed ASCI Supercomputer project,” SIGOPS Oper. Syst. Rev. 34(4), 76–96 (2000). [CrossRef]
45.	M. Pharr, C. Kolb, R. Gershbein, and P. Hanrahan, “Rendering Complex Scenes with Memory-Coherent Ray Tracing,” in SIGGRAPH ’97: Proceedings of the 24th annual conference on Computer graphics and interactive techniques (ACM Press/Addison-Wesley Publishing Co., 1997), pp. 101–108. [CrossRef]

OCIS Codes
(000.4430) General : Numerical approximation and analysis
(170.3660) Medical optics and biotechnology : Light propagation in tissues

ToC Category:
Medical Optics and Biotechnology

History
Original Manuscript: July 17, 2007
Revised Manuscript: September 2, 2007
Manuscript Accepted: September 3, 2007
Published: October 11, 2007

Virtual Issues
Vol. 2, Iss. 11 Virtual Journal for Biomedical Optics

Citation
Eduardo Margallo-Balbás and Patrick J. French, "Shape based Monte Carlo code for light transport in complex heterogeneous Tissues," Opt. Express 15, 14086-14098 (2007)
http://www.opticsinfobase.org/vjbo/abstract.cfm?URI=oe-15-21-14086

Sort: Author | Year | Journal | Reset

References

X. Gu, Q. Zhang, M. Bartlett, L. Schutz, L. Fajardo, and H. Jiang, "Differentiation of cysts from solid tumours in the breast with Diffuse Optical Tomography," Acad. Radiol. 11, 53-60 (2004). [CrossRef] [PubMed]
P. Taroni, G. Danesini, A. Torricelli, A. Pifferi, L. Spinelli and R. Cubbedu, "Clinical trial of time-resolved scanning optical mammography at 4 wavelengths between 683 and 975 nm," J. Biomed. Opt. 9, 464-473 (2004). [CrossRef] [PubMed]
S. Brown, E. Brown, and I. Walker, "The present and future role of photodynamic therapy in cancer treatment," Oncology 5, 497-508 (2004). [PubMed]
S. Chandrasekhar, Radiative Transfer (Oxford University Press, London, 1960).
L. Wang, S. Jacques, and L. Zheng, "MCML - Monte Carlo modeling of light transport in multi-layered tissues," Comput. Methods Programs Biomed. 47, 131-146 (1995). [CrossRef]
T. Pfefer, J. Barton, E. Chan, M. Ducros, B. Sorg, T. Milner, J. Nelson, and A. Welch, "A three-dimensional modular adaptable grid numerical model for light propagation during laser irradiation of skin tissue," IEEE J. Sel. Top. Quantum Electron. 2, 934-942 (1996). [CrossRef]
D. Boas, J. Culver, J. Stott, and A. Dunn, "Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head," Opt. Express 10, 159-170 (2002). [PubMed]
D. Côté and I. Vitkin, "Robust concentration determination of optically active molecules in turbid media with validated three-dimensional polarization sensitive Monte Carlo calculations," Opt. Express 13, 148-163 (2005). [CrossRef] [PubMed]
F. Brown, R. Barrett, T. Booth, J. Bull, L. Cox, R. Forster, J. Goorley, R. Mosteller, S. Post, R. Prael, E. Selcow, A. Sood and J. Sweezy, "MCNP Version 5," Trans. Am. Nucl. Soc. 87, 273 (2002).
J. Howell, "The Monte Carlo method in radiative heat transfer," J. Heat Transfer 120, 547-560 (1998). [CrossRef]
A. Marshak and A. B. Davis, eds., 3D Radiative Transfer in Cloudy Atmospheres (Springer-Verlag, Berlin, 2005). [CrossRef]
A. S. Glassner, Principles of Digital Image Synthesis (Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 1994).
E. Veach, "Robust Monte Carlo methods for light transport simulation," Ph.D. thesis, Stanford University (1997).
E. Cerezo, F. Pérez, X. Pueyo, F. Serón, and F. Sillion, "A survey on participating media rendering techniques," Visual. Comput. 21, 303-328 (2005). [CrossRef]
L. Szirmay-Kalos, V. Havran, B. Balász, and L. Szécsi, "On the efficiency of ray-shooting acceleration schemes," in SCCG ’02: Proceedings of the 18th spring conference on Computer graphics, pp. 97-106 (ACM Press, New York, NY, USA, 2002). [CrossRef]
A. Talsma, B. Chance, and R. Graaff, "Corrections for inhomogeneities in biological tissue caused by blood vessels," J. Opt. Soc. Am. A 18, 932-939 (2001). [CrossRef]
R. van Veen, W. Verkruysse, and H. Sterenborg, "Diffuse-reflectance spectroscopy from 500 to 1060 nm by correction for inhomogeneously distributed absorbers," Opt. Lett. 27, 246-248 (2002). [CrossRef]
J. Barton, T. Pfeffer, and A. Welch, "Optical Monte Carlo modeling of a true port wine stain anatomy," Opt. Express 2, 391-396 (1998). [CrossRef] [PubMed]
E. Okada, M. Firbank, M. Schweiger, S. Arridge, M. Cope, and D. 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]
Y. Fukui, Y. Ajichi, and E. Okada, "Monte Carlo prediction of near-infrared light propagation in realistic adult and neonatal head models," Appl. Opt. 42, 2881-2887 (2003). [CrossRef] [PubMed]
M. Firbank, S. Arridge, M. Schweiger, and D. Delpy, "An investigation of light transport through scattering bodies with non-scattering regions," Phys. Med. Biol. 41, 767-783 (1996). [CrossRef] [PubMed]
I. Patrikeev, Y. Petrov, I. Petrova, D. Prough, and R. Esenaliev, "Monte Carlo modeling of optoacoustic signals from human internal jugular veins," Appl. Opt., 46, 4820-4827 (2007). [CrossRef] [PubMed]
L. Murrer, H. Marijnissen, and W. Start, "Monte Carlo simulations for endobronchial photodynamic therapy: the influence of variations in optical and geometrical properties and of realistic and eccentric light sources," Laser Surg. Med. 22, 193-206 (1998). [CrossRef]
A. Pifferi, A. Torricelli, P. Taroni, A. Bassi, E. Chikoidze, E. Giambattistelli, and R. Cubeddu, "Optical biopsy of bone tissue: a step toward the diagnosis of bone pathologies," J. Biomed. Opt. 9, 474-479 (2004). [CrossRef] [PubMed]
M. Beers, R. Berkow, eds. The Merck Manual of Diagnosis and Therapy, 17th ed. (Merck Research Laboratories, Whitehouse Station, NJ, 1999).
F. Hartl, A. Tyndall, M. Kraenzlin, C. Bachmeier, C. Guckel, U. Senn, D. Hans, and R. Theiler, "Discriminatory ability of quantitative ultrasound parameters and bone mineral density in a population-based sample of postmenopausal women with vertebral fractures: results of the Basel Osteoporosis Study," J. Bone Miner. Res. 17, 321-330 (2002). [CrossRef] [PubMed]
S. Nayak, I. Olkin, H. Liu, M. Grabe, M. Gould, I. Allen, D. Owens, and D. Bravata, "Meta-Analysis: accuracy of quantitative ultrasound for identifying patients with osteoporosis," Ann. Intern. Med. 144, 832-841 (2006). [PubMed]
A. Takeuchi, R. Araki, S. Proskurin, Y. Takahashi, Y. Yamada, J. Ishii, S. Katayama, and A. Itabashi, "A new method of bone tissue measurement based upon light scattering," J. Bone Miner. Res. 12, 261-266 (1997). [CrossRef] [PubMed]
N. Ugryumova, S. Matcher, and D. Attenburrow, "Measurement of bone mineral density via light scattering," Phys. Med. Biol. 49, 469-483 (2004). [CrossRef] [PubMed]
J. Hebden, J. J. García Guerrero, V. Chernomordik, and A. H. Gandjbakhche, "Experimental evaluation of an anisotropic scattering model of a slab geometry," Opt. Lett. 29, 2518-2520 (2004). [CrossRef] [PubMed]
"HDF5 - A New Generation of HDF," Available at http://hdf.ncsa.uiuc.edu/HDF5 (2000). NCSA, University of Illinois at Urbana Champaign.
E. Margallo-Balbás, "Home page of TriMC3D," (2007), http://trimc3d.et.tudelft.nl/>
T. Möller and B. Trumbore, "Fast, minimum storage ray-triangle intersection," J. Graph. Tools 2, 21-28 (1997).
A. Glassner, "Space subdivision for fast ray tracing," IEEE Comput. Graphics Appl. 4, 15-22 (1984).
A. William, S. Barrus, R. Morley, and P. Shirley, "An Efficient and Robust Ray-Box Intersection Algorithm," J. Graph. Tools 10, 49-54 (2005). [CrossRef]
L. Wang, S. Jacques, and L. Zheng, "CONV - convolution for responses to a finite diameter photon beam incident on multi-layered tissues," Comput. Methods Programs Biomed. 54, 141-150 (1997). [CrossRef]
A. Tycho, T. Jørgensen, H. Yura, and P. Andersen, "Derivation of a Monte Carlo method for modeling heterodyne detection in optical Coherence tomography systems," Appl. Opt. 41, 6676-6691 (2002). [CrossRef] [PubMed]
E. Giesen and T. van Eijden, "The three-dimensional cancellous bone architecture of the human mandibular condyle," J. Dent. Res. 79, 957-963 (2000). [CrossRef] [PubMed]
M. Garland and P. Heckbert, "Surface simplification using Quadric Error Metrics," in SIGGRAPH (1997).
M. Firbank, M. Hiraoka, M. Essenpreis, and D. Delpy, "Measurement of the optical properties of the skull in the wavelength range 650-950 nm," Phys. Med. Biol. 38, 503-510 (1993). [CrossRef] [PubMed]
G. Bashkatov, "Optical properties of human skin, subcutaenous and mucous tissues in the wavelength range from 400 to 2000nm," J. Phys. D Appl. Phys. 38, 2543-2555 (2005). [CrossRef]
C. Tsai, Y. Yang, C. Han, J. Hsieh, and M. Chang, "Measurement and simulation of light distribution in biological tissues," Appl. Opt. 40, 5770-5777 (2001). [CrossRef]
R. van Veen, H. Sterenborg, A. Pifferi, A. Torricelli, E. Chikoidze, and R. Cubeddu, "Determination of visible near-IR absorption coefficients of mammalian fat using time and spatially resolved diffuse reflectance and transmission spectroscopy," J. Biomed. Opt. 10, 054004 (2005). [PubMed]
H. Bal, R. Bhoedjang, R. Hofman, C. Jacobs, T. Kielmann, J. Maassen, R. van Nieuwpoort, J. Romein, L. Renambot, T. Ruhl, R. Veldema, K. Verstoep, A. Baggio, G. Ballintijn, I. Kuz, G. Pierre, M. van Steen, A. Tanenbaum, G. Doornbos, D. Germans, H. Spoelder, E.-J. Baerends, S. van Gisbergen, H. Afsermanesh, D. van Albada, A. Belloum, D. Dubbeldam, Z. Hendrikse, B. Hertzberger, A. Hoekstra, K. Iskra, D. Kandhai, D. Koelma, F. van der Linden, B. Overeinder, P. Sloot, P. Spinnato, D. Epema, A. van Gemund, P. Jonker, A. Radulescu, C. van Reeuwijk, H. Sips, P. Knijnenburg, M. Lew, F. Sluiter, L. Wolters, H. Blom, C. de Laat, and A. van der Steen, "The distributed ASCI Supercomputer project," SIGOPS Oper. Syst. Rev. 34, 76-96 (2000). [CrossRef]
M. Pharr, C. Kolb, R. Gershbein and P. Hanrahan, "Rendering Complex Scenes with Memory-Coherent Ray Tracing," in SIGGRAPH ’97: Proceedings of the 24th annual conference on Computer graphics and interactive techniques (ACM Press/Addison-Wesley Publishing Co., 1997), pp. 101-108. [CrossRef]

Cited By

Alert me when this paper is cited

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

Click here to see a list of articles that cite this paper