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 , the scattering coefficient , and the anisotropy factor , which is the average value of the cosine of the scattering angles.4 A total interaction coefficient can also be defined. Its inverse is the mean free path between two consecutive interactions. Finally, if anisotropy is considered, the reduced scattering coefficient is defined asrepresents 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 , the time independent DA takes the form:being the average diffuse intensity, the source term, and 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 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. In that paper, the authors demonstrate that in the range of refractive indexes , 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 http://ifas.exa.unicen.edu.ar/es/grupos/j5-OpticaBiom/index.html.
The GPU we used is a GForce® 8600GT Nvidia (Santa Clara, California) graphic card installed in a PC with a Core 2 Duo, 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.
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.
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 -diam absorbing cylinder [ and ] immersed in a turbid slab with optical properties and . The cylinder extends vertically along the direction and is located on the optical axis with its longitudinal axis at 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.
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 ,6 will change with respect to the one when no inclusion is present , according tois the contribution of the inhomogeneity. With this notation, our proposal is to consider the ratio , resulting in the expression
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 -thick tank filled with a milky solution (see Sec. 3), the optical properties of which are given in Table 1 . A cylinder, high and 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.
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 1||Experiment 2||Experiment 3|
Please note that rotation of the optical axis (or the slab) by around the vertical axis produces a complementary view of the slab which, proceeding as described, gives the location of the inhomogeneity in the direction. For simplicity, in this example we have located the inhomogeneity at , 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 axis), and the comparisons are made along profiles taken at , 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.
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 and the reduced scattering coefficient simultaneously.5 In this technique, a pulsed laser working at repetition rate injects light pulses of very short duration (typically about ) 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:. 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 and -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 and an absorption coefficient , 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, laser diode operating at , and the images were collected by a , pixel, digital charge coupled device (CCD) camera and stored in the computer. A zoom lens was used to achieve a magnification ; thus, the resulting images map an area of the exit face of the tank. Figure 1 is a schema showing the experimental situation.
The initial MC simulations were carried out using 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 photons requires about 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 , giving an acceleration factor close to 600. However, after our first results were evaluated, it could be concluded that it is enough to use 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 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.
In experiment 1 (Figs. 4, 5, 6), the inclusion was a highly absorbing black painted wooden stick placed at different depths . 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.
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 1||Experiment 2||Experiment 3|
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.
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 -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 -long cylinder of diameter 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 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.
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.