Open Access
Add to BibSonomyAbstract
Conceptual engineering design and optimization of laser-based imaging techniques and optical diagnostic systems used in the field of biomedical optics requires a clear understanding of the light-tissue interaction and peculiarities of localization of the detected optical radiation within the medium. The description of photon migration within the turbid tissue-like media is based on the concept of radiative transfer that forms a basis of Monte Carlo (MC) modeling. An opportunity of direct simulation of influence of structural variations of biological tissues on the probing light makes MC a primary tool for biomedical optics and optical engineering. Due to the diversity of optical modalities utilizing different properties of light and mechanisms of light-tissue interactions a new MC code is typically required to be developed for the particular diagnostic application. In current paper introducing an object oriented concept of MC modeling and utilizing modern web applications we present the generalized online computational tool suitable for the major applications in biophotonics. The computation is supported by NVIDEA CUDA Graphics Processing Unit providing acceleration of modeling up to 340 times.
©2011 Optical Society of America
1. Introduction
Biophotonics and Biomedical Optics are well established and extremely important research areas traditionally associated with biomedical engineering and life sciences. A number of achievements in these fields arose due to a successful integration with optical diagnostic (OD) and optical/laser based imaging modalities. Nowadays OD technologies are widely used in a number of applications including tissue engineering, vascular and developmental biology, physiology, morphology, cancer research, material sciences, etc. [1–4]. OD provides a broad variety of practical solutions for clinical diagnostic and therapy in a range of studies from single cells and molecules to non-invasive biopsy of specific biological tissues and whole organs. The development of OD techniques for specific application is based on the practical implementation of laser/optical sources, detectors and various probe configurations. Conceptual design, further optimization and/or modification of the particular experimental OD systems requires careful selection of various technical parameters, including intensity, profile, coherence, polarization, wavelength of incident optical radiation; size, position, numerical aperture, sensitivity of the detector; as well as an estimation of how the spatial and temporal structural alterations in biological tissues can be linked with and can be distinguished by variations of these parameters. Therefore, an accurate realistic description of optical radiation propagation within complex biological media is required.
The description of optical radiation propagation within random media is based on the radiative transfer [5] that forms a basis of Monte Carlo (MC) modeling of photons migration in biological tissue [6]. Being originally introduced in biomedical optics for the counting of fluence rate distribution in biological tissues for the purpose of estimation laser radiation dose [7], in the last decades the MC approach has become a primary tool for a number of needs in biomedical optics. Incorporated with the computational model of human skin [8] MC technique has been used for simulation of skin optical & near-infrared reflectance spectra [9,10], analysis of skin fluorescence excitation [11,12], simulation of Optical Coherence Tomography (OCT) images of human skin [13,14], analysis of scattering orders and OCT image formation [15]. Using a combination of the MC technique and iterative procedure of the solution of Bethe-Salpeter equation, the MC approach has been generalized for simulation of coherent effects of multiple scattering, such as enhancement of Coherent Back-Scattering (CBS) and changes of temporal intensity correlation function depending on the dynamics of scattering particles [16,17]. Based on these developments a new approach of handling polarization has been introduced and some effects such as a helicity flip of circular polarization has been observed [18,19]. The obtained modeling results have been comprehensively validated by comparison with the known exact solution by Milne [20,21] and with the results of experimental studies of image transfer through the water solution of spherical micro-particles of known size and density [22,23]. Meanwhile, a number of other MC algorithms has been developed in the past, see for example [24–27].
It should be pointed out, however, that due to the diversity of optical modalities utilizing different properties of light and mechanisms of light-tissue interactions, including scattering, absorption, anisotropy, reflection, refraction, transmittance, interference, polarization, de-polarization, coherent scattering, phase shift, time-delaying, etc., a new MC code is typically required to be developed for the particular application. In current paper we introduce an object oriented concept of MC modeling designed as online computational tool for imitation of various aspects of optical/laser radiation propagation within complex disperse multi-layered media including biological tissues.
2. Object oriented concept of Monte Carlo modeling
Due to a number of practical OD applications the MC model undergoes continuous modifications and changes dedicated to inclusion of diverse properties of incident optical/laser radiation, configuration of the sources and detectors, structure of the medium and the conditions of light detection [8–27]. Past attempts to unify these MC codes are mainly based on the use of structured programming [28]. While structured programming is known for years, it limits the ability to handle a large code without decreasing its functionality and manageability [29]. In practice, the increasing diversity of the MC applications results to a substantial growth of the model’s source code and leads to the development of a set of separate MC codes dedicated each for a particular purpose.
To generalize and unify the code for a multi-purpose use in various biomedical optics applications we apply the Object Oriented Programming (OOP) concept. The OOP is widely used in mainstream application development and has been found extremely effective in design of complex multi-parametric systems, providing highly intuitive approach of programming [30]. The key features of OOP allow the MC to be separated into logical components, described by objects.
Thus, each photon packets has been defined as the object. The photon objects interact with the object - medium or medium components also defined as objects. Splitting the medium into the objects allows developing the tissues model more iteratively and uniformly. The distribution of scattering centers, macro-inhomogeneities, such as blood vessels, tumors, aneurisms, etc. can be formed by combination of 3D elementary volumes (objects) presenting spatial variations in the tissues. Moreover, actual structure of a biological tissue can be imported into the model as an object (image) provided by OCT, Photo-Acoustic Tomography (PAT), ultrasound, MRI, etc. The interaction between object - photons and object - medium originates similar as described in [6–23] and omitted for brevity in current paper.
Utilizing the inheritance feature of OOP the smart hierarchy structure of the code has been created to prevent creation of multiple classes for similar tasks. The hierarchy allows 'allied' objects to share variables and members, significantly reducing the amount of source code and paving the way to extend and generalize the MC for various OD applications. In addition, the variations of scattering phase functions, such as Mie, Rayleigh and Henyey-Greenstein [23] has been defined using the polymorphism feature of OOP that allows to handle the modeling with no changes of the source code.
Thus, the OOP approach significantly increases the efficiency of the model manageability and provides superior opportunities to generalize MC to combine previously developed MC models in a way to imitate a particular OD experiment taken into account various features of optical radiation and light-tissue interaction. Schematically the structure of the generalized MC approach is shown in Fig. 1 .
Fig. 1 Schematic presentation of the generalized object-oriented MC model structure. Data Input and selection of a particular application; Class objects represent exact experimental conditions including source-detector geometry & configuration, medium structure, properties of incident radiation, etc.
As one can see, first, by the selection of a certain application the parameters of the Source, Detector and Scattering medium are entered. Depending on the application, the objects are tuned to the appropriate feature of light-tissue interaction and the simulation is performed. By completing, the output data are prepared in the format utilized in the corresponding OD experiment/application. Current MC version supports the applications presented in [6–23] that includes estimation of sampling volume [8] offered by a reflectance and/or transmittance geometry (e.g. for fiber optic probe with a small ~mm source-detector spacing); simulation of reflectance or transmittance optical/near-infrared spectra of the multi-layered media like a human skin [9,10]; modeling of skin color depending on volume fraction of melanin and blood, blood oxygen saturation; modeling of OCT images with regards coherence and polarization properties of probing light (see the details in [13,14]), imitation of spatial localization of skin tissues fluorescence excitation [11,12] and simulation of simulation of coherent effects of multiple scattering [16,17].
Launching of a large number (~108) of photon packets and computing their interaction with medium and with the probe is a highly intensive computational process. Due to the task intensity, processing time has always been a significant issue in stochastic modeling, taking a few days to complete at the standard CPU. To achieve supreme performance of simulation a number programming approaches and optimizations of algorithms have been used in the past, including parallel and cluster computing [31,32].
We apply recently developed parallel computing framework, known as Compute Unified Device Architecture (CUDA), introduced by NVIDEA Corporation. NVIDEA CUDA technology provides an unlimited access to computational resources of graphic card: processor cores, different types of memory (of various capacity and speed) making Graphics Processing Unit (GPU) a massive co-processor in parallel computations [31,32].
We utilize the recent, introduced in 2010, CUDA generation, so-called architecture codename Fermi. Designed for C/C++ development and easily integrated with the Microsoft Visual Studio it makes parallel programming significantly easier, especially in terms of project management and debugging. The latest CUDA generation supports most of the OOP features like creating classes, inheritance, polymorphism and keeps all dramatic performance gains of the CUDA computing.
The OOP MC model has been developed using CUDA 3.2 C/C++ and supports multiple GPUs, providing user's interface and delivering the modeling output. The hardware is presented by two GeForce GTX 480 graphic cards with NVIDIA CUDA computing capability 2.0 totally having up to 2.7 Tflops of computational power on board (single precision). The graphic chip is logically divided in hundreds of CUDA cores (in our case 480 per card) schematically presented in Fig. 2 . Therefore, it executes up to twenty thousands of threads simultaneously, without context switch performance losses and very fast (up to 4 Gbit/s) on-chip GDDR memory. Graphic chip’s shared memory has been used to store the intermediate results; constant memory is applied for data input, whereas the global memory is used to store parameters of photon objects (e.g. path-length, state of polarization, outlet angles etc.). The tiling and cutoff techniques are used to process large data sets and avoid memory bandwidth limitations [31,32]. Such a design allows simulating propagation of thousands (15 to 20 depending on application and detector parameters) photon objects in the medium simultaneously that provides an opportunity for direct imitation of image or wave front transfer in the scattering medium [22] or OCT image modeling [13,14].
Fig. 2 Schematic presentation of used GPU logical divided on hundreds of independent cores allowing creation of thousands of lightweight parallel threads.
Specially designed for 3D graphic applications mainly for computer games and professional software designing this cutting edge graphic technology also incorporates a powerful set of instruments applied for optimized simulation of objects motion, rotation, reflection, ray-tracing, etc. The NVIDEA CUDA provides GPU-accelerated mathematical libraries, such as CULA, CUBLAS—Linear Algebra, CUFFT—Fast Furrier Transform, CURAND—Random Number Generators [33]. Their incorporation into MC allows speeding up the simulation of each photon packet up to 103 times.
3. Online solution and web integration
A conceptual design of the online solution is schematically presented on Fig. 3 .
Fig. 3 Schematic presentation of the key components of the online MC solution. The server hosts a web frontend, which accepts user simulation requests and retrieves results. The developed components provide interoperability between the interactive user interface and GPUs.
The MC modeling server provides major hardware acceleration by CUDA supporting GPUs. The developed server software consists of the web frontend, management part, GPU-Web integration and the CUDA-based MC core described above.
The frontend of the MC is developed by using Microsoft ASP.NET and Microsoft Silverlight technologies [34–37]. Microsoft Silverlight is used for creation of a cross-platform lightweight interactive interface to access a particular MC application [37], whereas ASP.NET is employed to meet modern design and security requirements [36].
The management part of the MC modeling server consists of the Input / Output (I/O) and Load Balancing systems. The I/O is created to accept and validate user credentials and online input of modeling parameters; to log and store computation enquires and output data into a server database; to deliver the results of simulation back to the frontend. The I/O works in a tandem with the Load Balancing which allows to manage the server load, monitor running simulations and to check the GPUs availability. To make the online MC less vulnerable to the common network threats and overloads, the management part has been developed utilizing recently added .NET 4.0 features, such as Windows Presentation Foundation (WPF) and Windows Communication Foundation (WCF) [38,39].
It should be pointed out that .NET solutions are executed in the Common Language Runtime (CLR), the core of .NET [35]. Known as a safe managed programming environment the CLR has no direct access to GPUs and other hardware resources. Therefore, to provide the interoperability between the web frontend and graphic cards the GPU-Web integration component has been developed. The Component Object Model (COM) [40] and .NET interoperability have been used to provide interaction between the binaries produced by C# and CUDA compilers. The proposed solution allows multiple users to access GPUs computational resources online within a reasonable time (in a range from a few seconds to several minutes).
The developed Online Object Oriented MC (O3MC) computational tool is now available at http://biophotonics.otago.ac.nz through the KAREN network [41], and can be used for various applications in biomedical optics as well as in biophotonics and optical engineering.
4. Results and discussion
Figure 4-a presents the interactive interface used for the selection of a particular application, and Figs. 4-b, 4-c, and 4-d show, respectively, the pages setting the parameters of the medium, probe configuration and the outline XYZ boundaries for sampling volume modeling.
Fig. 4 The frontend of O3MC computational tool: (a) interactive user interface providing the selection of applications; (b), (c) and (d) are, respectively, the web-pages that allow customizing parameters of the medium, fiber-optic probe and the boundaries of sampling volume. Available online at: http://biophotonics.otago.ac.nz/MCOnline.aspx.
The examples of results of sampling volume and skin reflectance spectra modeling counted by O3MC tool for a fiber-optic probe developed in house and used for the measurements of optical/near-infrared reflectance spectra of human skin in vivo are presented in Fig. 5 . Respectively, XZ (depth profile) and XY (surface) images of sampling volume for a particular position of source and detector are presented (see Fig. 5). The sampling volume presents a spatial localization of the detected signal within the skin tissues providing an estimation of which vascular bed is primarily responsible for the measured signal changes. The center-to-center source detector separation between fibers has been set as 400 μm, the beam profile has default as Gaussian that is typical for optical fibers. The source and detector areas can be overlapped or be coincided. The size and beam profile of the light source, numerical aperture, size and position of the detector are taken into account. YZ depth profile can be plotted as well.
Fig. 5 Example of sampling volume simulation for a fiber-optic probe applied for the measurements of human skin reflectance spectra by O3MC: side XZ (left) and surface XY (right) projections. Available online at http://biophotonics.otago.ac.nz/MCOnline.aspx.
Similar interface design and MATLAB style format of the results output are used in other applications included modeling human skin color, OCT images, transmittance spectra, CBS, confocal probing, etc. The results of simulation of spectra for multi-layered media and skin color are similar to the results presented in [9,10]; OCT images modeling and fluorescence are based, respectively, on [13,14] and [11,12]. Figure 6 presents an example of the results of human skin reflectance spectrum simulation and the results of skin color modeling with O3MC. The optical parameters of skin are taken same as reported in [9,10].
Fig. 6 Examples of reflectance spectrum simulation for Type I/II skin, melanin concentration 2% and normal concentration of blood (left) and skin color counting for Type V/VI skin, melanin concentration 25% and normal concentration of blood (right). Available online at http://biophotonics.otago.ac.nz/MCOnline.aspx.
The number of detected photon packets used in the simulation can be varied in a range of 105-108. Utilizing two GeForce GTX 480 graphic cards the O3MC requires 15.2 seconds to simulate the input of 108 photon packets. The time of simulation does not depend on absorption of the medium, but can be influenced by source-detector geometry increasing from a few seconds to 10 minutes. The results can be saved as an image or downloaded as a data file and then used in different scientific graphing software packages.
The developed O3MC approach has been comprehensively validated through the comparison of the results of diffuse reflectance Rd simulation with the known analytical results for semi-infinite medium tabulated by Giovanelli [42], with the results of adding-doubling method by Prahl [43] and with the results of MCML code developed by Wang et al. [7]. The following parameters of the medium have been used in the simulation: μs = 9 mm −1, μa = 1 mm −1, g = 0 (isotropic scattering), n = 1.5. The results of comparison are presented in Table 1 . Ten simulations with 106 photon packets each have been utilized in each case. Increasing number of photon packets to 108 reduces the error in O3MC for ~30%.
In addition, we compare the results of O3MC simulation of diffuse reflectance Rd and transmittance Td for a scattering slab illuminated by normally incident collimated light with the results tabulated by van de Hulst [44], with the results of adding-doubling method obtained by Prahl [43], and with the results of MCML developed by Wang et al. [7]. The following parameters of slab have been used in the simulation: μs = 9 mm −1, μa = 1 mm −1, g = 0.75, n = 1, slab thickness d = 0.2 mm. The results of comparison are listed in Table 2 .
The results presented in Tables 1 and 2 demonstrate a good accuracy of the developed code compare to known theoretical predictions and well agreed with the results of MC codes developed earlier [7,26,27]. With the further developments of O3MC the fluence rate, speckles, intensity/field correlation functions, fluorescence spectra, images of Photo-Acoustic Tomography and OCT, pulses propagation, image transfers and other applications currently used in biophotonics and biomedical optics will be added to online simulation at http://biophotonics.otago.ac.nz/MCOnline.aspx.
References and links
1. V. V. Tuchin, Handbook of Photonics for Biomedical Science (CRC Press, 2010).
2. D. A. Boas, C. Pitris, and N. Ramanujam, Handbook of Biomedical Optics (Taylor and Francis, 2010).
3. C. Balas, “Review of biomedical optical imaging—a powerful, non-invasive, non-ionizing technology for improving in vivo diagnosis,” Meas. Sci. Technol. 20(10), 104020 (2009). [CrossRef]
4. V. V. Tuchin, Handbook of Coherent Domain Optical Methods: Biomedical Diagnostics Environment and Material Science (Kluwer Academic, 2004).
5. L. S. Dolin, “Development of radiative transfer theory as applied to instrumental imaging in turbid media,” Phys.-Usp. 52(5), 519–526 (2009). [CrossRef]
6. F. Martelli, S. Del Bianco, A. Ismaelli, and G. Zaccanti, Light Propagation through Biological Tissue and Other Diffusive Media: Theory, Solutions, and Software (SPIE Press, 2009).
7. L. H. Wang, S. L. Jacques, and L.-Q. Zheng, “MCML—Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Methods Programs Biomed. 47(2), 131–146 (1995). [CrossRef]
8. I. V. Meglinsky and S. J. Matcher, “Modelling the sampling volume for skin blood oxygenation measurements,” Med. Biol. Eng. Comput. 39(1), 44–50 (2001). [CrossRef] [PubMed]
9. I. V. Meglinski, “Modeling the reflectance spectra of the optical radiation for random inhomogeneous multi-layered highly scattering and absorbing media by the Monte Carlo technique,” Quantum Electron. 31, 1101–1107 (2001).
10. I. V. Meglinski and S. J. Matcher, “Computer simulation of the skin reflectance spectra,” Comput. Methods Programs Biomed. 70(2), 179–186 (2003). [CrossRef] [PubMed]
11. D. Y. Churmakov, I. V. Meglinski, and D. A. Greenhalgh, “Amending of fluorescence sensor signal localization in human skin by matching of the refractive index,” J. Biomed. Opt. 9(2), 339–346 (2004). [CrossRef] [PubMed]
12. D. Y. Churmakov, I. V. Meglinski, S. A. Piletsky, and D. A. Greenhalgh, “Analysis of skin tissues spatial fluorescence distribution by the Monte Carlo simulation,” J. Phys. D Appl. Phys. 36(14), 1722–1728 (2003). [CrossRef]
13. I. V. Meglinski, M. Kirillin, V. L. Kuzmin, and R. Myllylä, “Simulation of polarization-sensitive optical coherence tomography images by a Monte Carlo method,” Opt. Lett. 33(14), 1581–1583 (2008). [CrossRef] [PubMed]
14. M. Kirillin, I. Meglinski, V. Kuzmin, E. Sergeeva, and R. Myllylä, “Simulation of optical coherence tomography images by Monte Carlo modeling based on polarization vector approach,” Opt. Express 18(21), 21714–21724 (2010). [CrossRef] [PubMed]
15. M. Y. Kirillin, I. V. Meglinskii, and A. V. Priezzhev, “Effect of photons of different scattering orders on the formation of a signal in optical low-coherence tomography of highly scattering media,” Quantum Electron. 36(3), 247–252 (2006). [CrossRef]
16. V. L. Kuzmin and I. V. Meglinski, “Coherent multiple scattering effects and Monte Carlo method,” JETP Lett. 79(3), 109–112 (2004). [CrossRef]
17. I. V. Meglinski, V. L. Kuzmin, D. Y. Churmakov, and D. A. Greenhalgh, “Monte Carlo simulation of coherent effects in multiple scattering,” Proc. R. Soc. A 461(2053), 43–53 (2005). [CrossRef]
18. V. L. Kuz’min and I. V. Meglinski, “Anomalous polarization effects during light scattering in random media,” J. Exp. Theor. Phys. 110(5), 742–753 (2010). [CrossRef]
19. V. L. Kuzmin and I. V. Meglinski, “Helicity flip of backscattered circularly polarized light,” Proc. SPIE 7573, 75730Z, 75730Z-7 (2010). [CrossRef]
20. V. L. Kuzmin and I. V. Meglinski, “Coherent effects of multiple scattering for scalar and electromagnetic fields: Monte-Carlo simulation and Milne-like solutions,” Opt. Commun. 273(2), 307–310 (2007). [CrossRef]
21. V. L. Kuz’min, I. V. Meglinski, and D. Y. Churmakov, “Stochastic modeling of coherent phenomena in strongly inhomogeneous media,” J. Exp. Theor. Phys. 101(1), 22–32 (2005). [CrossRef]
22. E. Berrocal, I. V. Meglinski, D. A. Greenhalgh, and M. A. Linne, “Image transfer through the complex scattering turbid media,” Laser Phys. Lett. 3(9), 464–467 (2006). [CrossRef]
23. E. Berrocal, D. Sedarsky, M. Paciaroni, I. V. Meglinski, and M. A. Linne, “Laser light scattering in turbid media Part I: Experimental and simulated results for the spatial intensity distribution,” Opt. Express 15(17), 10649–10665 (2007). [CrossRef] [PubMed]
24. 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 the adult human head,” Opt. Express 10(3), 159–170 (2002). [PubMed]
25. J. Ramella-Roman, S. Prahl, and S. Jacques, “Three Monte Carlo programs of polarized light transport into scattering media: part I,” Opt. Express 13(12), 4420–4438 (2005). [CrossRef] [PubMed]
26. Q. Fang, “Mesh-based Monte Carlo method using fast ray-tracing in Plücker coordinates,” Biomed. Opt. Express 1(1), 165–175 (2010). [CrossRef] [PubMed]
27. H. Shen and G. Wang, “A study on tetrahedron-based inhomogeneous Monte Carlo optical simulation,” Biomed. Opt. Express 2(1), 44–57 (2011). [CrossRef] [PubMed]
28. I. V. Meglinski, M. Kirillin, and V. L. Kuzmin, “The concept of a unified modelling of optical radiation propagation in complex turbid media,” Proc. SPIE 7142, 714204, 714204-5 (2008). [CrossRef]
29. S. Schach, Object-Oriented and Classical Software Engineering, 7th ed. (McGraw-Hill, 2006).
30. S. McConnell, Code Complete, 2nd ed. (Microsoft Press, 2004).
31. D. B. Kirk and W. W. Hwu, Programming Massively Parallel Processors: a hands-on Approach (MK Publishers, 2010).
32. J. Sanders and E. Kandrot, CUDA by Example: an Introduction to General-Purpose GPU Programming, (Addison-Wesley, 2010).
33. C. U. D. A. Programming Guide, 3.2. CUBLAS Library. CUFFT Library. CURAND Library, (NVIDIA Corporation, 2010).
34. L. Shklar and R. Rosen, Web Application Architecture: Principles, Protocols and Practices (Wiley, 2009).
35. A. Troelsen, Pro C# 2010 and the. NET 4 Platform (APRESS, 2010).
36. M. MacDonald, A. Freeman, and M. Szpuszta, Pro ASP.NET 4 in C# 2010, 4th ed. (APRESS, 2010).
37. L. Moroney, Microsoft Silverlight 4 Step by Step (Microsoft Press, 2010).
38. N. Pathak, Pro WCF 4: Practical Microsoft SOA Implementation (APRESS, 2011).
39. M. MacDonald, Pro WPF in C# 2010: Windows Presentation Foundation in. NET 4 (APRESS, 2010).
40. A. Troelsen, Developer's Workshop to COM and ATL 3.0 (Wordware Publishing, 2000).
41. See details at “Kiwi advanced research and education network,” http://karen.net.nz/home/.
42. R. Giovanelli, “Reflection by semi-infinite diffusers,” Opt. Acta (Lond.) 2, 153–162 (1955).
43. S. Prahl, “The adding-doubling method,” in Optical-Thermal Response of Laser Irradiate Tissue, A. Welch and M. van Gemert, eds. (Plenum, 1995), Chap.5, pp. 101–125.
44. H. van de Hulst, Light Scattering by Small Particles (Dover, 1957).
