Translator Disclaimer
1 May 2010 Solution of the direct problem in turbid media with inclusions using Monte Carlo simulations implemented in graphics processing units: new criterion for processing transmittance data
Author Affiliations +
The study of light propagation in diffusive media requires solving the radiative transfer equation, or eventually, the diffusion approximation. Except for some cases involving simple geometries, the problem with immersed inclusions has not been solved. Also, Monte Carlo (MC) calculations have become a gold standard for simulating photon migration in turbid media, although they have the drawback large processing times. The purpose of this work is two-fold: first, we introduce a new processing criterion to retrieve information about the location and shape of absorbing inclusions based on normalization to the background intensity, when no inhomogeneities are present. Second, we demonstrate the feasibility of including inhomogeneities in MC simulations implemented in graphics processing units, achieving large acceleration factors (~103), thus providing an important tool for iteratively solving the forward problem to retrieve the optical properties of the inclusion. Results using a cw source are compared with MC outcomes showing very good agreement.



The study of near-infrared (NIR) light propagation in turbid media has been of great interest for scientists during the last decade because of several potential applications, with optical imaging in biological tissue being of utmost importance. Using this noninvasive radiation has led to the development of many biomedical applications, from detection of breast cancer1 to bedside monitoring of cerebral blood flow,2, 3 to mention but a few.

Propagation of light inside turbid media can be completely described in terms of several parameters, namely the absorption coefficient μa , the scattering coefficient μs , and the anisotropy factor g=cosθ , which is the average value of the cosine of the scattering angles.4 A total interaction coefficient μt=μs+μa can also be defined. Its inverse μt1 is the mean free path between two consecutive interactions. Finally, if anisotropy is considered, the reduced scattering coefficient is defined as

where (μs)1 represents the transport mean free path.4

From a theoretical point of view, light propagation in highly scattering media is described by the scalar radiative transfer equation (RTE),4 which takes into account the energy balance inside a volume of the scattering medium. However, since the RTE is relatively complicated, in many practical cases dealing with media where scattering dominates over absorption, the diffusion approximation (DA)4 is preferred. In particular, for an isotropic source and for μsμa , the time independent DA takes the form:

Eq. 1

with U(r) being the average diffuse intensity, E(r) the source term, and D=13(μs+μa) is the diffusion coefficient.

Even though analytical solutions of the DA for homogeneous media can be obtained with relative ease,5, 6 many biomedical applications of interest, such as optical mammography, require solving the problem considering immersed inhomogeneities; the solution becomes more complicated in this case and it is still matter of active research.7

Monte Carlo (MC) simulations, in which photons are tracked while they propagate inside the medium, have been used as a gold standard to validate simple models in diffusion problems; it is a very flexible, powerful tool that allows one to obtain a solution for almost any imaginable geometry and configuration, and it is thus capable of dealing with relatively complex situations.

In general, the desired solution is the diffusively transmitted or reflected number of photons. Thus, one needs to accumulate at the output those photons that have succeeded in propagating through the turbid medium and results are built up stochastically. Depending on the optical properties of the medium, it is possible that a large number of photons must be launched before one is collected, a condition that normally requires very long calculation times to achieve results with reasonable good statistics.8, 9, 10, 11, 12 This individual tracking of photons, despite the fact that the recent development of faster and/or multiple-core processors has outstandingly reduced calculation time, remains without doubt the main drawback of MC simulations, and has motivated scientists to seek for alternative acceleration procedures that range from tracking only useful (or successful photons)13 to implementing hybrid algorithm mixing theory with pure numerical simulations.14 However, a real breakthrough in reducing time in MC simulations can be found in very recent publications involving the parallel processing capabilities of graphics processor units (GPUs). 15, 16, 17

On the other hand, despite the importance of taking into account immersed inhomogeneities in turbid media, there have been only a few contributions dealing with numerical simulations considering diffusive inclusions inside diffusive media.14, 18

In this work we compare MC simulations to CW light transmittance profiles of turbid media slabs with inclusions. However, these profiles are not always clear indicators of which kind of inhomogeneity is immersed, nor where it is placed. To overcome this difficulty, we propose a new processing criterion based on normalization relative to the transmittance profiles of the homogeneous medium without inclusions. As it is shown, this approach allows retrieving the location of the inhomogeneity as well as its approximate shape and dimensions. With this previous information available, the MC code is run for several sets of optical parameters until a good match to the experimental profiles is found, thus resulting in a “most likely” set of optical parameters for the inclusion.

To overcome the drawback of the intrinsic time-consuming nature of MC simulations, we used a very fast MC algorithm for simulations in turbid media containing inhomogeneities that performs computing on GPU instead of using the CPU, and is based on compute unified device architecture (CUDA). Due to the parallel processing capability of GPU, and depending on the graphics card used, calculation times can be drastically reduced. This is a very important feature, since, even though the essential inverse problem of finding the inclusion that produces a given light distribution at the output cannot be addressed, the previous information about location of the inclusion and the great flexibility of MC, combined with the high speed of parallel processing, give the possibility of multiple trials in reasonable time. Moreover, when testing any inversion algorithm, it is desirable to have a set of simulated “tomographic images” for which the exact location of the inhomogeneity as well as its actual optical parameters are precisely known. These images can be obtained in a relative short time using parallel processing in the MC calculations.

In particular, we have modified the open source code CUDAMCML for layered turbid media presented in Ref. 15, and which can be found in Ref. 19 by giving the possibility of including inhomogeneities. The multiple layer structure proposed in Ref. 19 is kept and inhomogeneities (one or more) can be located inside the desired layer. Our modified code treats the case of inhomogeneities embedded only in one layer. Photons are launched normal to the entrance face of a slab of thickness S and they are tracked until they are either absorbed, or exit the medium at the opposite side within a defined detection area. Fresnel reflections at both flat slab boundaries, namely entrance and exit face, are considered as in the original code. However, taking into account Fresnel reflections at the bulk - inclusion boundary can be very time consuming, since when a photon reaches this interface, the incidence angle must be evaluated, which in turn depends on the chosen geometry for the inclusion. Instead, in our proposal, these reflections are considered according to the suggestion of Ripoll and Nieto Vesperinas,20 in which the actual reduced scattering coefficient of the inclusion is replaced by an effective one defined by

being nrel=nbulkninclusion . In that paper, the authors demonstrate that in the range of refractive indexes 1.3<n<1.5 , as it is typical for biological tissues, the error committed in the worse case is never above 3%.

The code resulting from this contribution can be freely accessed at

The GPU we used is a GForce® 8600GT Nvidia (Santa Clara, California) graphic card installed in a PC with a Core 2 Duo, 2.66-GHz CPU. Even for this card, which is clearly not a last generation one, we measured acceleration factors in the order of several hundreds.

In short, this contribution is organized as follows. In the next section we give a detailed presentation containing the conceptual basis of the proposed normalization criterion. Next, in Sec. 3, we describe and compare experiments with the corresponding MC simulations. Main results are also discussed in this section. Finally, the most important conclusions of the work are presented.


Analysis Criterion

As shortly mentioned in Sec. 1, we present a criterion for retrieving additional information about the inclusions, which, to the best of our knowledge, has not been reported yet. Based on theoretical solutions for simple cases such as spheres21 and cylinders22, 23 embedded in slabs, we found that the ratio between light transmittance when the inhomogeneity is present to that without inclusion provides a good analysis criterion.

In fact, as it is shown in the experiments, the transmittance profiles do not always allow one to deduce the presence of an inclusion. In contrast, the proposed analysis criterion is capable of retrieving the optical properties of the inclusion. It is also sensitive to the location of the inhomogeneity, which constitutes valuable a-priori information to be used in the MC simulations.


Conceptual Basis of the Proposed Approach

When using CW light, two different situations are possible to explore a slab containing a turbid medium with inclusions.

  • 1. Transilluminance is when the laser and detector face each other (defining an optical axis), and are moved together scanning the slab along a line (1-D case) or plane (2-D case). Detection is made over a small area at the exit face of the slab by an optical fiber placed in contact with it or by imaging a small area with a camera. Thus, the detected light has explored a relative small volume inside the slab. This volume contains the most probable photon paths inside the slab, and defines the so-called photon hitting density, a figure that is also often known as the “banana.” The presence of an inclusion produces variations in the intensity at some points of the scan. A plot of intensity versus position gives rise to a figure which is self “normalized” with respect to the baseline intensity (a position for the laser-detector pair, where the banana does not touch the inhomogeneity). Lateral resolution is dependent on both the slab thickness and the optical properties of the bulk, since they affect the shape of the “banana” inside the slab. This procedure can give good lateral resolution but, as a drawback, it requires, as mentioned, scanning the laser and detector along a line or a matrix of points, with the consequent time consumption.

  • 2. Transmittance is the case presented in this work accordingly with the setup of Fig. 1 . No scanning is made and a single picture of an extended area at the light output is acquired. Using the concept of photon hitting density, it can be said that the intensity at each point of this picture is the result of the interaction of multiple “bananas” defined between the source and the multiple points at the image. Data acquisition is now much easier, since no translation of the laser-detector pair is required (no scanning). However, as a drawback, the resulting picture is not normalized and it does not always provide evident information about the presence or absence of an inclusion.

Fig. 1

Sketch of the experimental setup (top view). A CW laser diode, operating at 830nm , illuminates one face of the tank containing the milky diffusive solution and the inclusion. The collimated laser beam defines the optical axis of the system, which is coincident with the optical axis of the camera lens. A positioning device (not shown) allows placing the inclusion in the transverse direction (x,y) and at different depths along the optical axis (z) .


The proposed criterion is the result of fusing together these two approaches: a single picture is acquired that is simple and fast, and after that it is normalized to the background transmittance without inclusion, as happens in transilluminance.

As a proof of principles, we show in Figs. 2, 2, 2 the results obtained when using the theoretical model of a cylinder immersed in a slab of turbid media22, 23 for the two cases mentioned before. For this example we present in Fig. 2 the transilluminance and transmittance profiles, with no normalization, for a 1.4-cm -diam absorbing cylinder [ μa(cyl)=0.28cm1 and μs(cyl)=10cm1 ] immersed in a 3-cm turbid slab with optical properties μa=0.08cm1 and μs=10cm1 . The cylinder extends vertically along the y direction and is located on the optical axis (x=0) with its longitudinal axis at z=1.5cm depth (see Fig. 1 for the geometry). Note that the relative higher absorption of the cylinder is not obvious in the transmittance curve. Figure 2 shows the transmittance profile for the slab with the cylindrical inclusion of Fig. 2, together with the transmittance of the slab without inclusion (background transmittance). Figure 2 compares the quotient of the two curves of Fig. 2, which is the ansatz of this work, with the transilluminance profile of Fig. 2.

Fig. 2

Conceptual basis of this proposal. (a) Theoretical transmittance and transilluminance profiles when a cylindrical inclusion is present inside a slab. See Sec. 2.1 for the corresponding optical parameters used. (b) Theoretical transmittances with inclusion [J(r)] and without inclusion [J0(r)] ; note that J(r) alone does not gives a clear indication about the optical properties of the inhomogeneity. (c) Normalized fluxes: transilluminance normalized to the baseline value and transmittance normalized as proposed in this work [see Eq. 3].


As can be seen from Fig. 2, both curves tend to the same value far from the optical axis; therefore, the background transmittance becomes a natural “normalizer” for the transmittance profile. In fact, this quotient clearly shows the presence of the more absorbing inclusion, and this ratio has the same qualitative behavior as the transilluminance profile. For comparison, this last curve has been normalized to its own background (far from the optical axis). Thus, both constructions show the presence of the absorbing inclusion. Note that the normalized transilluminance profile has a better lateral resolution, while the modulation is similar to that of the normalized transmittance. However, as was already said, the construction of a transilluminance profile is much more time consuming, since it requires scanning the laser-detector pair.

From a formal point of view, because of the introduction of the inclusion, the photon flux at the detector position, which using conventional notation is given by J(r)=4πDϕ(r) ,6 will change with respect to the one when no inclusion is present J0(r) , according to

Eq. 2

where Jinh(r) is the contribution of the inhomogeneity. With this notation, our proposal is to consider the ratio Q(r)=J(r)J0(r) , resulting in the expression

Eq. 3

therefore accounting for the presence of the inhomogeneity.


Real Situations


Optical properties of the bulk

It is very important to remark that the application of the proposed approach requires knowing the transmittance of the bulk without inclusions (background transmittance). This is an obvious drawback if trying to apply it to real situations, where inclusions cannot be removed.

There are several ways of dealing with this situation. 1. The first rough approach is to use the optical properties published for several kinds of typical healthy tissues3 and, given a slab thickness, to obtain, via MC simulations or theoretical models for homogeneous media slabs,6 the corresponding background transmittance.

2. A better solution to this problem consists in measuring the optical properties at several different positions of the slab for which, in principle, it is not known if it contains inclusions or not (the detailed procedure for measuring optical properties is discussed later in Sec. 3.1). Thus, provided that the inclusion, if present, represents a perturbation (as is usually the case), the average values of the optical properties measured at these points will give representative values of the optical properties of the bulk without inclusions. Even though this solution may differ from the actual case considered, using this background would provide a starting situation giving a first approximation at least about the location of the inhomogeneity (inhomogeneities). This information can, in a next step, be used in the MC simulations to proceed to the retrieval of the optical properties of the inclusion. This is shown later in the experiments.


Location of the inclusion

Before MC simulations are carried out to retrieve the optical properties of the inclusion, its location (in 3-D) must be inferred. As an example of how this works in our proposal, we present in Fig. 3 an experiment for a S=3-cm -thick tank filled with a milky solution (see Sec. 3), the optical properties of which are given in Table 1 . A cylinder, 16mm high and 16mm in diameter, is placed inside the tank as an absorbing inclusion. Figure 3 is the complete transmittance image, with the slab containing the inclusion. A background transmittance image is also acquired without the inclusion. Please refer to Sec. 3 for the details about image acquisition. Figure 3 is the result of normalizing the transmittance image of Fig. 3 to the actual experimental transmittance background. Note that the location of the inhomogeneity (unknown a priori) relative to the optical axis results in Fig. 3. The horizontal and vertical profiles shown help to better visualize the location of the inclusion. In Fig. 3 the experimental background (which in real situations is not available) is replaced by one generated by MC using average optical properties of the slab. This is applied to all experiments presented in Sec. 3.2. Despite the fact that the modulations in the profiles of Fig. 3 are not the same as those in Fig. 3, the primary information concerning the location of the inclusion and its “more absorbing” nature are identical in both figures.

Fig. 3

2-D images, as seen by the camera, and the corresponding horizontal and vertical profiles to show how our proposal works. (a) A complete transmittance image of the tank containing an absorbing inclusion with unknown location; profiles are taken at the optical axis ( x=0 and y=0 ). (b) The result of normalizing the image shown in (a) to the actual experimental background (without inclusion); this was done in this controlled experiment to compare with the result shown in (c). (c) Same as (b) but normalized to a background generated by MC simulations using the optical properties given in Table 1. Please note that in (b) and (c), the location of the center of the inclusion is practically the same, found to be at x18mm and y2mm , relative to the optical axis (0,0). Additionally, a clue about the lateral extension of the inclusion can be inferred from the FWHM of these profiles.


Table 1

Measured optical properties of the bulk and inclusions (within 10%) corresponding to the different experiments. All measurements were made following the procedure described in Sec. 3.1.

Experiment 1Experiment 2Experiment 3
μa[cm1] 0.07 0.07
μs[cm1] 9.049.409.409.409.40

Please note that rotation of the optical axis (or the slab) by 90deg around the vertical axis produces a complementary view of the slab which, proceeding as described, gives the location of the inhomogeneity in the z direction. For simplicity, in this example we have located the inhomogeneity at z=S2 , and thus only one view of the tank is needed to have all three coordinates of the center of the inclusion.


Optical properties of the inclusion

Once the location of the inclusion is known, we propose to exploit the high speed calculation capabilities of MC implemented in GPU to vary the optical properties of the inclusion to find the best match to the experiment. In accordance to this proposal, we present in the next section a set of controlled experiments that illustrate the capability of retrieving optical properties of inclusions by selecting the MC simulation that produces the best match. Since direct comparison of 2-D images can be relatively cumbersome, all experiments presented refer to cylindrical inhomogeneities (extended along the y axis), and the comparisons are made along profiles taken at y=0 , which can be directly compared in a similar figure. All transmittance profiles are normalized to a background generated by MC using the optical properties of the bulk listed in Table 1.


Experimental Methods


Measurement of Optical Properties

In the experiments presented in this work, measurement of optical properties is mandatory. To this end, we used a time-resolved approach as described in the following. The advantage of this time-resolved method is that it is possible to retrieve both the absorption coefficient μa and the reduced scattering coefficient μs simultaneously.5 In this technique, a pulsed laser working at 50-MHz repetition rate injects light pulses of very short duration (typically about 70ps ) in an optically turbid medium via an optical fiber. Photons are diffused in the medium, collected at the exit by another optical fiber, and sent to a high sensitivity detector [photomultiplier tube (PMT)]. Due to dispersion inside the medium, photon trajectories are different from each other, and thus they last different times to reach the collecting optical fiber; this produces a temporal distribution of photons depending on the optical properties of the medium, which is called distribution of times of flight (DTOF). It is built up and displayed on a monitor by means of specific electronic equipment, which is the basis of the technique known as time correlated single photon counting (TCSPC), in which photons are assumed to arrive one by one at the detector. A detailed description of this technique can be found in the monograph by Becker.24 Fitting the resulting DTOF with an adequate theoretical model allows us to obtain the desired optical properties; since we are dealing with slabs of turbid media, we use the model described by Contini, Martelli, and Zaccanti.6

It must be mentioned that due to the transit time of photons along the optical fibers and inside the detector, even if no turbid medium is present the laser pulses are temporally spread, giving rise to a certain instrument response function (IRF), which is also measured experimentally by placing the laser fiber face to face with that going to the PMT. The measured DTOF is then the convolution of this IRF with the actual desired temporal distribution; this must be taken into account when the fitting procedure is carried out. In symbols:

where the desired function is DTOFactual . Since deconvolution can be numerically unstable, the usual procedure is to convolve the measured IRF with a theoretical DTOF6, 25 for a set of optical parameters and to compare this result with the measured DTOF. The process is iterated in a Levenberg−Marquardt algorithm for different sets of optical properties until the desired convergence is achieved.


Experiments, Simulations, and Discussion

To illustrate how our proposal works, we present three controlled experiments for CW light transmittance in turbid media slabs, in which diffusive light intensity, experimentally measured or simulated, is mapped on a bidimensional array at the exit face of the slab. From these 2-D intensity arrays (images), 1-D profiles are taken at a constant vertical position for better visualization and simpler processing. Inclusions consist of cylindrical inhomogeneities placed at different depths inside a homogeneous slab. Optical properties of both, the bulk and the inhomogeneities for all three experiments, are summarized in Table 1. They were measured using the method described in the previous section. This constitutes the a-priori information to test the goodness of the approach.

Experiments were carried out using a 25×25-cm and 3-cm -thick cuvette filled with a solution of milk, distilled water, and 1:1000 diluted india ink in proportions 1:3:0.04. This mix gives a reduced scattering coefficient (μs) and an absorption coefficient (μa) , which are very close to those of biological tissues. They were measured previous to the transmittance experiments by fitting DTOF of photons obtained by TCSPC, as described before. The inclusions, when present, were held in place inside the milky solution by a rigid mount attached to a translation stage for proper positioning. Illumination was provided by the collimated beam of a CW, 10-mW laser diode operating at λ=830nm , and the images were collected by a 14-bit , 768×512 (horizontal×vertical) pixel, digital charge coupled device (CCD) camera and stored in the computer. A zoom lens was used to achieve a magnification m=0.1 ; thus, the resulting images map an area 76.8×51.2mm2 of the exit face of the tank. Figure 1 is a schema showing the experimental situation.

The initial MC simulations were carried out using 1010 photons. This is important in transmittance experiments; at the wings of the profiles, far from the optical axis, the number of photons is poor and decreases quickly. Simulating the paths of 1010 photons requires about 30to40min using the GPU. This time must be compared to the one taken by the Core 2 Duo CPU for performing the same number of simulations, which is about 12days (1.7×104min) , giving an acceleration factor close to 600. However, after our first results were evaluated, it could be concluded that it is enough to use 109 photons in each simulation; thus, calculation time is reduced to a few minutes. So, it is possible to obtain a representative set of MC outcomes (say ten) in a reasonable lapse, which allows us to infer the optical parameters of the inclusion, as well as re-estimate its size and location. Accordingly to this, all MC results presented in Figs. 4, 5, 6, 7, 8, 9, 10 correspond to simulations with 109 photons. Clearly the exact number of photons required varies accordingly with both the optical properties and the thickness of the slab and the type of inclusion.

Fig. 4

Highly absorbing inclusion (black painted wooden stick), ϕinc=0.7cm , placed next to the light entrance face of the tank. (a) Normalized experimental intensity profile (circles) traced at vertical pixel 256 (y=0) , shown together with the normalized outcome of the MC simulation (solid curve). (b) Experiment and MC profiles given in (a) divided by a common background generated by MC simulations, accordingly with the proposal explained in the text.


Fig. 5

Same as in Fig. 4 for the highly absorbing inclusion placed at the middle of the tank.


Fig. 6

Same as in Fig. 4 for the highly absorbing inclusion placed next to the light exit face of the cuvette.


Fig. 7

Absorbing agarose cylinder of ϕinc=1.4cm placed next to the light entrance face of the tank. (a) Normalized experimental intensity profile (circles) traced at vertical pixel 256 (y=0) , shown together with the normalized outcome of the MC simulation (solid curve). (b) Experiment and MC profiles given in (a) divided by a common background generated by MC simulations, accordingly with the proposal explained in the text.


Fig. 8

Same as in Fig. 7 for the agarose cylinder at the middle of the tank.


Fig. 9

Same as in Fig. 7 for the agarose cylinder next to the light exit face of the tank.


Fig. 10

Small diameter agarose cylinder less absorbing than the surrounding medium placed next to the light exit face of the tank. (a) Normalized experimental intensity profile (circles) traced at vertical pixel 256 (y=0) , shown together with the normalized outcome of the MC simulation (solid curve). (b) Experiment and MC profiles given in (a) divided by a common background generated by MC simulations, accordingly with the proposal explained in the text.


In experiment 1 (Figs. 4, 5, 6), the inclusion was a highly absorbing black painted wooden stick (ϕinc=0.7cm) placed at different depths z . In all cases, part (a) refers to the normalized transmittance profiles with inclusion obtained by both experiment and MC simulation. Part (b) shows the results after applying the analysis of Sec. 2.2.2. As stated, profiles in (b) are obtained after division by a background generated by MC, since in general, the experimental background is not known.

The important features that simulations are expected to reproduce are both full width at half maximum (FWHM) and intensity modulation of the intensity profiles. That means that those profiles obtained from MC must match the corresponding experimental ones. When this condition is reached, the best set of optical parameters of the inclusion has been found. These retrieved values for the experiments presented next are summarized in Table 2 . Note that intensity modulation as well as FWHM of the transmittance profile are both expected to change for different depths of the same inclusion. This is because the most probable paths for the photons (the “banana”) have no constant section inside the tank, and thus the number of photons intersected by the inclusion changes for different depths of it.

Table 2

Retrieved optical properties of the inclusions corresponding to the different experiments according to the proposal of this work. The given uncertainties are not of statistical nature; they arise because for each experiment, three positions of the inclusion must be considered simultaneously.

Experiment 1Experiment 2Experiment 3
μa[cm1] > 10 0.23±0.05 0.020±0.002
μs[cm1] 9.40±1.00 9.40±1.00

In Fig. 4, the inclusion was placed inside the tank next to the face where light enters the slab. Both FWHM and intensity modulation retrieved from MC agree with the experimental data, as can be directly seen from Fig. 4. Notice that this figure is ambiguous about the optical properties of the inclusion with respect to the bulk. However, after processing [Fig. 4], it becomes evident that an inclusion, which is more absorbing than the bulk, is present.

Figures 5 and 6 are analogous to Fig. 4 but with the highly absorbing inclusion placed in the middle of the tank (z=S2) and close to the wall facing the camera, respectively.

Please note the ambiguity of Figs. 5 and 6, which at first sight could lead to the conclusion that two inclusions that are less absorbing than the bulk are present inside the tank. This ambiguity is solved in both cases after normalization to the background, as shown in Figs. 5 and 6. Note also that, as expected, Figs. 4, 5, 6 show that the more distant the inclusion from the camera, the smaller the intensity modulation and the wider the profile.

As another test much closer to an actual situation, we present experiment 2, where the inclusion is a cylinder made of agarose mixed with a milky solution similar to that of the tank and india ink to achieve the desired absorption. The optical properties of the inclusion could not be measured directly from the inclusion itself because of its relatively small dimensions and its shape. Instead, we used the same mixture of agarose and milky solution to fabricate, simultaneously with then inclusion, a 3-cm -thick slab from which we measured the optical properties by fitting the DTOF of photons, as described before. Since the process of solidification of agarose inside the small plastic cylindrical form used to fabricate the inclusion may differ from the corresponding one when generating the slab (much greater volume), slight departures of the actual optical properties of inclusion from those of the slab could be expected. Results for this case are presented in Figs. 7, 8, 9. For this series of experiments, we used the same tank filled with the same milky solution as experiment 1. The inclusion, a 25-cm -long cylinder of diameter ϕinc=1.4cm had an absorption coefficient four times higher than that of the solution filling the tank. Refractive index of both bulk and inclusion were assumed equal.

Again, the comparisons between experiment and simulations for transmittance data and processed profiles show very good agreement. Sensitivity with depth is also demonstrated from the series of Figs. 7, 8, 9, which show changes in the slenderness of the profiles.

As a final example, we present in experiment 3 (Fig. 10) the case of a slightly perturbing less absorbing inclusion, consisting of an agarose cylinder of ϕinc=0.5cm immersed in the cuvette filled with the same milky solution used for the preceding experiments. It was placed next to the light exit face of the slab. Notice that the processed profile [Fig. 10] produces a peak that is higher than the background, clearly indicating the presence of a less absorbing inclusion.

Note that in Figs. 4, 5, 6, 7, 8, 9, 10, all transverse positions are measured with respect to the optical axis defined by the line joining the laser and the camera (Fig. 1).



The present work deals with the retrieval of (previously known) optical properties by iteratively using high speed MC simulations. Two general conclusions can be mentioned. 1. It demonstrates that normalization of a transmittance image to a transmittance background using MC gives the same information as a transilluminance experiment, which is much more complicated, allowing one to obtain both the position and approximate dimensions of immersed inhomogeneities. A clue about the optical properties of the inclusion relative to the bulk (more or less absorbing) is also obtained, thus simplifying the next step. 2. Given the results of 1., the high calculation speed of MC implemented in GPU permits one to retrieve the optical properties of the inclusion by iteratively solving the forward problem.

In particular, it is demonstrated from the analysis summarized in Figs. 2 and 3 that simple observation of a transmittance image alone does not give a clear indication about the relative absorption of an inclusion. However, after dividing by the transmittance background, the relative absorption of the inclusion becomes evident, as is the case in transilluminance experiments. Since transmittance is performed in 2-D, and transilluminance requires scanning the surface of the medium under study, our proposal greatly simplifies the data acquisition due to the consequent time reduction.

Our proposal is based on the equivalence of Figs. 3 and 3 when the normalization transmittance background is obtained by an experiment without inclusion [Fig. 3] or by MC [Fig. 3], the same essential information is obtained. This is very important for real cases where the actual experimental image without inclusions cannot be obtained.

The method unambiguously indicates if the inclusion is more or less absorbing than the bulk, even for the case of slightly perturbing conditions. Moreover, in three series of controlled experiments (Figs. 4, 5, 6, 7, 8, 9, 10), we show that the optical properties of the inclusion are retrieved with good accuracy, as is shown in Table 2.


The authors would like to acknowledge the financial support from grant PICT 38058. Carbone acknowledges the financial support from a CICPBA scholarship. We would also like to thank H. F. Ranea Sandoval for useful suggestions and critical reading of the manuscript.



A. Pifferi, P. Taroni, A. Toricelli, F. Messina, R. Cubeddu, and G. Danesini, “Four wavelength time resolved optical mammography in the 680–980 range,” Opt. Lett., 28 (13), 1138 –1140 (2003). 0146-9592 Google Scholar


A. Liebert, H. Wabnitz, J. Steinbrink, H. Obrig, M. Möller, R. Macdonald, A. Villringer, and H. Rinneberg, “Time-resolved multidistance near-infrared spectroscopy of the adult head: intracerebral and extracerebral absorption changes from moments of distribution of times of flight of photons,” Appl. Opt., 43 (15), 3037 –3047 (2004). 0003-6935 Google Scholar


Biomedical Photonics Handbook, CRC Press, Boca Raton, FL (2003). Google Scholar


A. Ishimaru, Wave Propagation and Scattering in Random Media, IEEE Press, Oxford University Press, Oxford, UK (1997). Google Scholar


M. S. Patterson, B. Chance, and B. C. Wilson, “Time resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties,” Appl. Opt., 28 (12), 2331 –2336 (1989). 0003-6935 Google Scholar


D. Contini, F. Martelli, and G. Zaccanti, “Photon migration through a turbid slab described by a model based on diffusion approximation. I. Theory,” Appl. Opt., 36 (19), 4587 –4599 (1997). 0003-6935 Google Scholar


S. D. Konecky, G. Y. Panasyuk, K. Lee, V. Markel, A. G. Yodh, and J. C. Schotland, “Imaging complex structures with diffuse light,” Opt. Express, 16 (7), 5048 –5060 (2008). 1094-4087 Google Scholar


D. A. Boas, J. P. Culver, J. J. Stott, and A. K. Dunn, “Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including adult human head,” Opt. Express, 10 (3), 159 –170 (2002). 1094-4087 Google Scholar


G. S. Fishman, Monte Carlo Concepts, Algorithms and Applications, Springer, New York (1996). Google Scholar


F. Martelli, D. Contini, A. Taddeucci, and G. Zaccanti, “Photon migration through a turbid slab described by a model based on diffusion approximation. II. Comparison with Monte Carlo results,” Appl. Opt., 36 (19), 4600 –4612 (1997). 0003-6935 Google Scholar


H. Key, E. R. Davies, P. C. Jackson, and P. N. T. Wells, “Monte Carlo modelling of light propagation in breast tissue,” Phys. Med. Biol., 36 (5), 591 –602 (1991). 0031-9155 Google Scholar


E. B. de Haller and C. Depeursinge, “Simulation of time-resolved breast transillumination,” Med. Biol. Eng. Comput., 31 (2), 165 –170 (1993). 0140-0118 Google Scholar


A. Sassaroli, C. Blumetti, F. Martelli, L. Alianelli, D. Contini, A. Ismaelli, and G. Zaccanti, “Monte Carlo procedure for investigating light propagation and imaging of highly scattering media,” Appl. Opt., 37 (31), 7392 –7400 (1998). 0003-6935 Google Scholar


H. O. Di Rocco, D. I. Iriarte, J. A. Pomarico, and H. F. Ranea-Sandoval, “Acceleration of Monte Carlo modeling of light transport in turbid media; an approach based on hybrid, theoretical and numerical calculations,” J. Quant. Spectrosc. Radiat. Transf., 110 (4–5), 307 –314 (2009). 0022-4073 Google Scholar


E. Alerstam, T. Svensson, and S. Andersson-Engels, “Parallel computing with graphics processing units for high-speed Monte Carlo simulation of photon migration,” J. Biomed. Opt., 13 (6), 060504 (2008). 1083-3668 Google Scholar


Q. Fang and D. A. Boas, “Monte Carlo simulation of photon migration in 3-D turbid media accelerated by graphics processing units,” Opt. Express, 17 (22), 20178 –20190 (2009). 1094-4087 Google Scholar


W. C. Lo, K. Redmond, J. Luu, P. Chow, J. Rose, and L. Lilge, “Hardware acceleration of a Monte Carlo simulation for photodynamic treatment planning,” J. Biomed. Opt., 14 (1), 014019 –014030 (2009). 1083-3668 Google Scholar


L. Dagdug, V. Chernomordik, G. H. Weiss, and A. H. Grandjbakhche, “Monte Carlo simulations of increased/decreased scattering inclusions inside a turbid slab,” Phys. Med. Biol., 50 (23), 5573 –5581 (2005). 0031-9155 Google Scholar


J. Ripoll and M. Nieto-Vesperinas, “Index mismatch for diffuse photon density waves at both flat and rough diffuse-diffuse interfaces,” J. Opt. Soc. Am. A, 16 (8), 1947 –1957 (1999). 0740-3232 Google Scholar


D. A. Boas, M. A. O’Leary, B. Chance, and A. G. Yodh, “Scattering of diffuse photon density waves by spherical inhomogeneities within turbid media: analytic solution and applications,” Proc. Natl. Acad. Sci. U.S.A., 91 4887 –4891 (1994). 0027-8424 Google Scholar


S. A. Walker, D. A. Boas, and E. Gratton, “Photon density waves scattered from cylindrical inhomogeneities: theory and experiments,” Appl. Opt., 37 (10), 1935 –1944 (1998). 0003-6935 Google Scholar


H. O. Di Rocco, D. I. Iriarte, J. A. Pomarico, and H. F. Ranea-Sandoval, “CW laser transilluminance in turbid media with cylindrical inclusions,” Optik, Google Scholar


W. Becker, Advanced Time-Correlated Single Photon Counting Techniques, Springer Verlag, Berlin (2005). Google Scholar


S. J. Matcher, “Closed form expresions for obtaining the absorption and scattering coefficients of a turbid medium with time-resolved spectroscopy,” Appl. Opt., 36 (31), 8298 –8302 (1997). 0003-6935 Google Scholar
©(2010) Society of Photo-Optical Instrumentation Engineers (SPIE)
Nicolas Carbone, Hector Di Rocco, Daniela I. Iriarte, and Juan A. Pomarico "Solution of the direct problem in turbid media with inclusions using Monte Carlo simulations implemented in graphics processing units: new criterion for processing transmittance data," Journal of Biomedical Optics 15(3), 035002 (1 May 2010).
Published: 1 May 2010

Back to Top