1 March 2013 Monte Carlo lookup table-based inverse model for extracting optical properties from tissue-simulating phantoms using diffuse reflectance spectroscopy
Author Affiliations +
We present a Monte Carlo lookup table (MCLUT)-based inverse model for extracting optical properties from tissue-simulating phantoms. This model is valid for close source-detector separation and highly absorbing tissues. The MCLUT is based entirely on Monte Carlo simulation, which was implemented using a graphics processing unit. We used tissue-simulating phantoms to determine the accuracy of the MCLUT inverse model. Our results show strong agreement between extracted and expected optical properties, with errors rate of 1.74% for extracted reduced scattering values, 0.74% for extracted absorption values, and 2.42% for extracted hemoglobin concentration values.



Diffuse reflectance spectroscopy (DRS) has been widely used to characterize tissue optical properties for disease diagnosis.12.3.4 Typically, DRS uses a fiber to deliver light to tissue. The delivered light is both scattered and absorbed by the tissue and is then recollected by another fiber a short distance from the source fiber. The collected light, or diffuse reflectance, contains quantitative information about tissue structure and composition. While the instrumentation for a DRS system is very simple, the accurate extraction of optical properties from the collected signal is a significant challenge. An accurate tissue model is needed to relate the collected signal to tissue optical properties. One method for analyzing diffuse reflectance spectra relies on the solution to the diffuse approximation of the radiative transport equation. However, the diffusion approximation is not valid for short source-detector separations or in highly absorbing tissues.5

Because most cancers originate in the epithelial layer, it is necessary to use probes with a short source-detector distance in order to sample photons that travel through the epithelial layer.6 In addition, angiogenesis, an indicator of early cancer, can cause a significant increase in absorption due to blood. Unfortunately, the diffuse approximation is not valid in these regimes. Recently, reflectance lookup tables have been used to analyze diffuse reflectance spectra.7,8 These lookup tables are created in two different ways: experimental measurements of phantoms with known optical properties7 or Monte Carlo simulations.9 Creating a lookup table (LUT) with experimental measurements has the advantage of incorporating unknown system responses into the LUT. However, the creation of an experimental LUT is time-consuming, and the accuracy is dependent on the skill and experience of the investigator. The Monte Carlo method is especially useful for creating an LUT, because it provides the ability to model complex probe geometries and tissue structures. However, intensive computation is required to achieve results with desirable variance, which can make it extremely time-consuming to populate an LUT containing thousands of values.

Much prior work has been undertaken to improve the speed and efficiency of the Monte Carlo method for modeling light transport in turbid media. These methods can be separated into three groups:

  • 1. methods that use the information from a small set of Monte Carlo simulations and scale the results to calculate a wide range of optical properties10

  • 2. methods that use geometry splitting techniques to increase the fraction of useful photons11

  • 3. methods that parallelize the Monte Carlo simulations.12

The first set of methods has the advantage of not requiring a large number of simulations to create an LUT. However, errors arise because the photon trajectory information necessary to perform the scaling operation can be recorded only in several depth intervals with finite widths.13 Geometry splitting techniques decrease the number of photons needed for a simulation to converge, but their implementation is very difficult for complex probe geometries. The third set of methods uses classical Monte Carlo simulations and therefore does not make any sacrifices in accuracy, flexibility, or implementation difficulty. Speedup is achieved by simulating multiple photons simultaneously on different processors. Because each photon is independent of every other photon, this problem is considered “embarrassingly parallel.” Alerstam et al. have shown that general-purpose graphic processing units (GPGPUs) can increase the speed of Monte Carlo simulations of photon transport by three orders of magnitude on a relatively inexpensive GPU when compared to the sequential implementation.12 We present, for the first time, a Monte Carlo LUT-based model where all values in the LUT were created by independent Monte Carlo simulations by using a parallel implementation on a GPGPU.


Creation of the Lookup Table

A two-dimensional Monte Carlo code written in ANSI C implemented on an NVIDIA GTX 560 Ti GPU with NVIDIA’s Compute Unified Device Architecture (CUDA) was used to simulate photon reflectance in a single-layer tissue model on 386 parallel threads.12 The multiply with carry random number generator was used. The refractive index above the tissue was set to 1.452 to match the refractive index of the fiber, and the refractive index of the tissue was set to 1.33 to match the refractive index of the phantoms. To test the effect of errors in the refractive index, we created LUTs with different refractive indices and repeated the extraction of optical properties with the different LUTs. We found a 5% error rate in the refractive index corresponds to error increases of 1.3% for extracted μs values and 0.8% for extracted μa values. To prevent photons from exiting the tissue volume, the radius and width of the tissue volume were set to 3 cm. A total of 1×106 photons were launched to obtain the impulse response. To ensure stochastic noise would be sufficiently low in LUT locations with high albedo, 100 separate MC simulations were performed with the optical properties in the LUT that would give the lowest value of reflectance. We found that using 1×106 photons reduced the standard deviation of the 100 different reflectance values to less than 0.5% of the mean. The diffuse reflectance for our specific probe geometry was then calculated by convolving the impulse response with the beam profile.14 Our probe was modeled using a Gaussian shaped beam profile of collimated light with a diameter of 200 μm, a detector diameter of 200 μm, and a source-detector separation of 250 μm. The diffuse reflectance values for all physiological realistic combinations of scattering and absorption were calculated using the GPGPU Monte Carlo implementation with the tissue and probe geometry described above. Twenty increments (N) were used for both scattering (0 to 50cm1) and absorption (1 to 50cm1), meaning a total of 400 separate Monte Carlo simulations were needed to create the LUT. The Henyey-Greenstein phase function was used for sampling scattering angles. The scattering anisotropy (g) was set to 0.85 for all simulations. For the range of g values present in human tissue (g>0.8), it has been shown10 that the diffuse reflectance will be the same for any values of μsand g that generate the same μs. The resulting MCLUT is shown in Fig. 1. It took 2 min to run the 400 separate Monte Carlo simulations.

Fig. 1

The resulting lookup table [R(μs,μa)] created using 400 separate Monte Carlo simulations. Each Monte Carlo simulation was used to calculate a reflectance value for a given scattering coefficient and absorption coefficient combination.



Forward and Inverse Models

For the forward model used to generate diffuse reflectance spectra, the reduced scattering coefficient was constrained to the form


where λ0=630nm. The absorption coefficient was calculated using


where εi(λ) is the wavelength-dependent extinction coefficient of a chromophore, Ci is the concentration of that chromophore, and N is the number of chromophores. Depending on the type of tissue sampled and the wavelength range of interest, any number of chromophores can be used to calculate μa(λ). Once μs(λ) and μa(λ) are calculated, the MCLUT can be used to generate a modeled reflectance spectrum. Cubic splines were used to interpolate between values in the LUT. Figure 2(a) shows the forward model of diffuse reflectance.

Fig. 2

Flowcharts for the forward (a) and inverse (b) models of diffuse reflectance used to create the modeled spectrum and to fit the MCLUT model to the reflectance data.


Figure 2(b) shows the inverse model used to fit our diffuse reflectance spectra. First, we made an initial guess for the optical properties, and then the forward model was used to generate a spectrum. Next, the sum of squares error between the predicted reflectance and the measured reflectance was calculated using


where δ is the sum of squares error, K is the number of wavelength points, Rm is the measured spectrum, Re is the modeled spectrum, and λi is the wavelength. The parameters are then iteratively updated until the sum of squares error is minimized. An interior-point nonlinear optimization routine provided in the MATLAB optimization toolbox (Mathworks, Nattick, Massachusetts) was used as the optimization algorithm. The average fit time was 1.3 s, with the number of iterations limited to 5000 and the termination tolerance on the error function set to 1×103. Because the modeled spectra are in absolute units (photons counted) and the measured spectra are in units relative to a baseline calibration measurement, it is necessary to perform a calibration so that the modeled and measured spectra can be compared. Additionally, the calibration corrects for wavelength-dependent responses in the experiment that are not accounted for in our forward model. We performed the calibration by taking the ratio of a modeled spectrum and a measured spectrum with the same optical properties. Then, to make the measured spectra equivalent to the modeled spectra, all measured spectra were multiplied by this ratio. To ensure that the choice of optical properties used in the calibration step did not bias the results, the mean of the ratios for all spectra used in the validation set was used for calibration. The mean calibration ratio is shown in Fig. 3(a).

Fig. 3

(a) A diffuse reflectance spectrum (dashed) [μs(λ0)=25.4cm1 and [Hb]=1.5mg/ml] and the associated MCLUT-fit (solid). (b) Mean of all ratios of measured and modeled diffuse reflectance spectra with the same optical properties. This ratio was used for calibration.



Validation and Results

To test the performance of our MCLUT-based inverse model, we created 21 tissue phantoms with hemoglobin (Hb) (Sigma-Aldrich) as the absorber and polystyrene beads (diameter=1μm) as the scatterer. Hb concentration [(Hb)] ranged from 0 to 3mg/ml, and the reduced scattering coefficient [μs(λ0=630nm)] ranged from 6.4 to 27.5cm1. We used Mie theory to calculate the μs' of the tissue phantoms. For our inverse model, we assumed the absorption in the visible range was due to oxy-hemoglobin. We measured the optical density for the HbO2 solution using a spectrophotometer and calculated the absorption spectrum using Beer’s law. Because the addition of HbO2 dilutes the solution, a small change in μs was accounted for when calculating the known values for μs. The DRS system consisted of a xenon flash lamp (Model: E6611, Hamamatsu) as the light source, a spectrograph (Model SP2150i, Princeton Intruments) and camera (Cool-SNAP, Photometrics) as the spectrometer; and a fiber optic probe with the geometry described above (FiberTechOptica, Ontario, Canada). The diffuse reflectance spectrum and its associated fit can be seen in Fig. 3(b) and shows that the inverse model can accurately fit the experimental data. Figure 4(a) and 4(b) shows scatter plots of the extracted versus expected for μs(λ0) and [Hb], respectively. The solid line in each graph is the line of perfect agreement. The results indicate there is excellent agreement between the extracted and expected optical properties. The MCLUT inverse model estimated the optical properties over a wide range with average root-mean-square percent errors of 1.74% for μs, 0.74% for μa, and 2.42% for [Hb]. We compared the performance of our MCLUT-based model to an experimental LUT-based model. The MCLUT model was able to estimate μs and μa with decreases in percent error of 3.16% and 10.86%, respectively, when compared to the experimental LUT model.7

Fig. 4

(a) Hemoglobin concentration extracted from the MCLUT inverse model versus known hemoglobin concentration. The solid line indicates perfect agreement. (b) μs(λ0) extracted from the MCLUT inverse model versus known μs(λ0). The solid line indicates perfect agreement. (c) In vivo reflectance spectra from two representative groups: clinically normal and BCC. The thin sold line is the fit.


Our inverse model was then tested on previously collected data from a clinical feasibility study15 to illustrate the application of the model for noninvasive detection of skin cancer. Figure 4(c) shows representative spectra from two groups: clinically normal and basal cell carcinoma (BCC). The plot shows good agreement between the MCLUT fit and the measured in vivo spectra. For this analysis, the absorption coefficients calculated in the forward model were determined by melanin, deoxy-hemoglobin, and oxy-hemoglobin.



Although other inverse models have recently been developed for extracting optical properties from diffuse reflectance spectra, this work represents the first Monte Carlo lookup table-based inverse model where the LUT was generated entirely by Monte Carlo simulation. Previously, the amount of time required to generate a LUT entirely by Monte Carlo simulation made this technique infeasible. However, recent advances in GPGPU computing have allowed parallel Monte Carlo implementations capable of running three orders of magnitude faster than traditional, serial implementations of Monte Carlo simulations. By creating an LUT entirely by Monte Carlo simulation, our method is not subject to the errors that arise from using either the diffusion approximation or the Monte Carlo scaling method. When compared to an experimental LUT, our method was more accurate, but, more importantly, it has the advantages of being easier to implement and repeatable. This model can also be adapted to more complex probe and tissue geometries.


The funding for this project was provided by the NIH (NIBIB R21EB0115892 and NCI RO1 CA132032), the NSF (DGE-1110007), and the Texas Higher Education Coordinating Board.


1. G. ZoniosJ. BykowskiN. Kollias, “Skin melanin, hemoglobin, and light scattering properties can be quantitatively assessed in vivo using diffuse reflectance spectroscopy,” J. Invest. Dermatol. 117(6), 1452–1457 (2001).JIDEAE0022-202X http://dx.doi.org/10.1046/j.0022-202x.2001.01577.x Google Scholar

2. G. Zonioset al., “Diffuse reflectance spectroscopy of human adenomatous colon polyps in vivo,” Appl. Opt. 38(31), 6628–6637 (1999).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.38.006628 Google Scholar

3. B. W. Murphyet al., “Toward the discrimination of early melanoma from common and dysplastic nevus using fiber optic diffuse reflectance spectroscopy,” J. Biomed. Opt. 10(6), 064020 (2005).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.2135799 Google Scholar

4. U. UtzingerR. Richards-Kortum, “Fiber optic probes for biomedical optical spectroscopy,” J. Biomed. Opt. 8(1), 121–147 (2003).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.1528207 Google Scholar

5. T. FarrellM. PattersonB. Wilson, “A diffusion-theory model of spatially resolved, steady-state diffuse reflectance for the noninvasive determination of tissue optical-properties invivo,” Med. Phys. 19(4), 879–888 (1992).MPHYA60094-2405 http://dx.doi.org/10.1118/1.596777 Google Scholar

6. Y. N. Mirabalet al., “Reflectance spectroscopy for in vivo detection of cervical precancer,” J. Biomed. Opt. 7(4), 587–594 (2002).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.1502675 Google Scholar

7. N. RajaramT. H. NguyenJ. W. Tunnell, “Lookup table-based inverse model for determining optical properties of turbid media,” J. Biomed. Opt. 13(5), 050501 (2008).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.2981797 Google Scholar

8. B. NicholsN. RajaramJ. Tunnell, “Performance of a lookup table-based approach for measuring tissue optical properties with diffuse optical spectroscopy,” J. Biomed. Opt. 17(5), 057001 (2012).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.JBO.17.5.057001 Google Scholar

9. G. M. PalmerN. Ramanujam, “Monte Carlo-based inverse model for calculating tissue optical properties. Part I: Theory and validation on synthetic phantoms,” Appl. Opt. 45(5), 1062–1071 (2006).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.45.001062 Google Scholar

10. R. Graaffet al., “Condensed Monte Carlo simulations for the description of light transport,” Appl. Opt. 32(4), 426–434 (1993).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.32.000426 Google Scholar

11. E. TinetS. AvrillierJ. Tualle, “Fast Semianalytical Monte Carlo simulation for time-resolved light propagation in turbid media,” JOSA A 13(9), 1903–1915 (1996).JOAOD60740-3232 http://dx.doi.org/10.1364/JOSAA.13.001903 Google Scholar

12. E. AlerstamT. SvenssonS. Andersson-Engels, “Parallel computing with graphics processing units for high-speed Monte Carlo simulation of photon migration,” J. Biomed. Opt. 13(6), 060504 (2008).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.3041496 Google Scholar

13. Q. LiuN. Ramanujam, “Scaling method for fast Monte Carlo simulation of diffuse reflectance spectra from multilayered turbid media,” JOSA A 24(4), 1011–1025 (2007).JOAOD60740-3232 http://dx.doi.org/10.1364/JOSAA.24.001011 Google Scholar

14. L. WangS. JacquesL. Zheng, “CONV—convolution for responses to a finite diameter photon beam incident on multi-layered tissues,” Comput. Methods Programs Biomed. 54(3), 141–150 (1997).CMPBEK0169-2607 http://dx.doi.org/10.1016/S0169-2607(97)00021-7 Google Scholar

15. N. Rajaramet al., “Pilot clinical study for quantitative spectral diagnosis of non-melanoma skin cancer,” Lasers Surg. Med. 42(10), 716–727 (2010).LSMEDI0196-8092 http://dx.doi.org/10.1002/lsm.21009 Google Scholar

© The Authors. Published by SPIE under a Creative Commons Attribution 3.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.
Ricky J. Hennessy, Ricky J. Hennessy, Sam L. Lim, Sam L. Lim, Mia K. Markey, Mia K. Markey, James W. Tunnell, James W. Tunnell, "Monte Carlo lookup table-based inverse model for extracting optical properties from tissue-simulating phantoms using diffuse reflectance spectroscopy," Journal of Biomedical Optics 18(3), 037003 (1 March 2013). https://doi.org/10.1117/1.JBO.18.3.037003 . Submission:

Back to Top