Monte Carlo method is applied for simulation of 2D optical coherence tomography (OCT) images of skin-like model. Layer boundaries in skin model feature curved shape which agrees with physiological structure of human skin. The effect of coherence properties of probing radiation on OCT image formation and speckles in the detected OCT signal is considered. The developed model is employed for image simulation both for conventional and polarization dependent time-domain OCT modalities. Simulation of polarized OCT signal is performed using vector approach developed previously for modeling of electromagnetic field transfer in turbid media.
©2010 Optical Society of America
Nowadays optical coherent tomography (OCT) is well known as a reputable diagnostic technique for non-invasive imaging of complex turbid media such as biological tissues [1–5]. Nevertheless, a wide application of OCT in day-to-day clinical practice is significantly limited due to image blurring and in-depth signal decrease caused by multiple scattering of probing radiation within the biological tissues . To improve the OCT system design and image processing algorithms a systematic understanding of OCT image formation is required.
Due to a greater complication arisen when polarization and coherence properties of scattered laser radiation should be taken into account the radiation transfer cannot be described simply in terms of diffusion theory. Moreover, for a complex multi-layered media with non-planar layer boundaries encountered in practical applications it is rarely possible to find an analytical solution. This situation is typical for many practical applications including OCT. Numerical techniques such as Monte Carlo (MC) method are typically applied. MC technique is well established in various biomedical [7,8], astronomical and meteorological applications [9–11], industrial sprays  and others. A number of studies has been devoted, in the past, to analyze contribution of scattering orders into the image transfer [12,13] and formation of OCT signal [8,14–19].
In current report, we introduce MC-based technique applied for simulation of 2D OCT images of skin-like model with curved layer boundaries corresponding to the actual skin structure confirmed by OCT images of normal human skin obtained in vivo.
2. Experimental OCT imaging of skin
The principles of OCT are widely described elsewhere [1–5]. Experimental study was performed using the OCT system developed at the Institute of Applied Physics RAS (Nizhny Novgorod, Russia) [20,21]. OCT setup operates at the central wavelength of 910 nm with spectrum width of 50 nm corresponding to axial and transversal spatial resolutions of 15 and 25 μm, respectively. The examples of typical OCT images of human skin obtained in vivo by conventional (Fig. 1a ) and polarization dependent (Fig. 1b and Fig. 1c) modalities exhibit complex multi-layered tissues structure.
These images clearly show high scattering in upper Stratum Corneum layer manifested by bright area in non- and co-polarized modes. However, its brightness in the cross-polarization image is significantly smaller because of its small thickness which does not provide sufficient depolarization. Next layer, lower Stratum Corneum, seems relatively dark in all three images, whereas underlying layer epidermis appears bright in the obtained OCT images in non- and co-polarized modes (see Fig. 1a and 1b, correspondingly). High scattering in this layer also results in significant depolarization degree which is evident by high intensity of this layer in the cross-polarization image (see Fig. 1c).
3. Numerical simulation
3.1 Basic concept of Monte Carlo simulation of OCT signal
Principles of stochastic MC method for numerical calculation of radiation intensity scattered within a randomly inhomogeneous turbid medium are widely described in the literature [7–9,12–19]. MC is based on the consequent simulation of a random photon trajectory within the medium between the point where the photon enters the medium, and the point where it leaves the medium. Simulation of the photon trajectories consists of the following key stages: injection of the photon in the medium, generation of the photon path-length, generation of a scattering event, definition of reflection/refraction at the medium boundaries, definition of detection and accounting for the absorption. The photon free path s between the two successive elastic scattering events is determined by the Poisson probability density function :
The probability that the photon free path exceeds s is defined as:
Given the probability density function (1) it is easy to express the random magnitude s via the probability ξ:
Direction of the photon after single scattering is defined by the scattering phase function:23].
Described steps are repeated till the photon is detected arriving at the detector with the given area and acceptance angle, or leaves the scattering medium, or its total path exceeds the maximum allowed path. The total number of the launched photons used in the simulation is typically ~107. The details of the reflection and refraction at the medium boundary and at the interface between layers are given in .26]:18]:
In order to simulate the 2D OCT image consequent A-scans are simulated with the definite transversal step in probing position. The total number of simulated A-scans and the transversal step between them are predefined regarding the width (FWHM) of the probing beam diameter.
3.2 Simulation of polarization dependent OCT signal
Polarization of an electromagnetic wave is typically described in framework of the Stokes-Mueller or Jones formalism [27,28]. Stokes-Mueller formalism was applied to study polarization in birefringent turbid media for potential application to polarization-sensitive optical imaging . The experimental and numerical studies were carried out to observe the backscattering polarization patterns presented in a form of the Mueller matrices [30,31]. The residual polarization degree of the backscattered light and its connection to the optical properties of the scattering medium was studied in . To explore the possibility of retrieving the birefringence properties of layered tissue with the depth-resolved polarization-sensitive OCT (PS-OCT) Jones formalism was implemented to reduce large calculation time, typically required in the full Stokes-Mueller approach .
Schmitt et al. proposed a method for discriminating short and long path photons that is based on the relationship between the polarization states of incident and forward scattered light . Discussing coherent backscattering and its dependence on the state of polarization of an incident linearly polarized light, Akkermans et al. assumed that for both states of polarization the corresponding intensity can be represented as a product of the intensity of a scalar wave, which does not depend on the polarization state, and a corresponding multiplicative factor (weighting function) describing polarization transfer . These factors are specific for a given polarization state and generally depend on the properties of scattering particles, e.g. size, shape [36,37]. The expression for the co- and cross-polarized intensities derived from the expression proposed in  has been successfully used in , whereas the experimental results  proved the adequacy of the Akkermans’ conjecture of time scales involved into depolarization process of the backscattered light.
Consider a plane electromagnetic wave polarized in the x direction that enters the medium along positive direction of z axis normal to the interface. By a co-polarized wave we understand a linearly polarized scattered wave with the same orientation of polarization as the incident wave, and a cross-polarized wave is perpendicular to the incident wave direction of polarization . Thus, waves scattered in xz and yz planes define co-polarized and cross-polarized components of the scattered electromagnetic wave, respectively.
To account for depolarization effect in PS-OCT images we adopted the polarization vector formalism where the polarization is described in terms of a polarization vector undergoing a sequence of transformations after each scattering event . The trajectories of the polarized photons are weighted in accordance with their polarization state. Within the far-field or Fraunhofer approximation the polarization vector of the scattered wave is transformed upon the i-th scattering event into as [23,35]:23]. Namely, the size D of the particles should obey the relation (εr-1)D/λ << 1 where (εr-1) is the relative fluctuation of dielectric permittivity between the scatterer (e.g. cell component such as nucleus or mitochondria) and the surrounding medium (e.g. cytoplasm). Typically in biotissues the value of (εr −1) is less than 0.1  therefore RGD approximation is quite reasonable for the particles with the sizes of units of λ which are characterized by non-isotropic scattering phase function. Recently it has been demonstrated that the expression (10) can be successfully applied with the use of Heyney-Greenstein phase function . At the same time, implementation of Eq. (10) for strongly scattering large inclusions may lead to some discrepancy in the final results; however, this approximation can be applied to qualitatively estimate the effects of depolarization in OCT images. Explicitly, the tensor is presented as38,42,43]:
Finally, for linearly polarized probing radiation the expression (9) can be re-written separately for co- and cross-polarized detected OCT signal:
4. Modeling of OCT images of skin
Basing on the principles of cross-polarization OCT signals formation described above we performed MC simulation of the OCT images of skin in non-, co- and cross-polarized modes corresponding to those presented in Fig. 1. Instrumental parameters used in simulation are taken in accordance with the characteristics of the reference experimental cross-polarization OCT setup described in section 2. The incidence of radiation at the surface of turbid medium is considered normal. This situation corresponds to the case of collimated probing beam. The case when imaging is performed with the focused Gaussian-shaped probing beam the waist of which is localized inside the medium is typically simulated using-ray-optics approach proposed by Milsom . For simulation of an A-scan 5·106 photons are employed; 50 A-scans are obtained to construct each 2D OCT image.
A sophisticated non-planar multilayered model of skin based on a number of experimental OCT images of skin of a male thumb was used in simulations. It is typically characterized by thick Stratum Corneum layer which can be separated into two parts: upper thin layer consisting of chaotically oriented dry cells and lower thick layer consisting of ordered cells (Fig. 2 ). The underlying skin layers are represented by epidermis, dermis with upper plexus, dermis and dermis with lower plexus . The shape of the boundaries and the thickness of the layers are chosen basing on the observed experimental OCT images. The optical properties of the layers are taken from  with a slight modification giving the larger weight to the values extracted from the experimental OCT images which mostly suits the simulation conditions. The chosen optical properties represented in Table 1 cannot be considered as exact ones due to significant variations of reported values from subject to subject and depending on extraction method, thus they can serve only as approximate values.
5. Results and discussion
Simulated OCT images of skin in non-polarized mode implemented using the expression (9) without account for “speckle effects” are presented in Fig. 3 for four various values of coherence lengths. The coherence length of the probing radiation varies from l c = 5 up to 30 μm at the wavelength λ = 910 nm that is typical for super-luminescence diodes employed in OCT setups. All figures are shown in the same color scale, so the increase in coherence length results in increase of corresponding average signal intensity. At the same time the axial resolution of the images decreases with the increase of coherence length thus yielding blurring of the image. The obtained images qualitatively depict the essential features of the experimental images (Fig. 1a) exhibiting bright epidermis layer below dark lower Stratum Corneum layer. Such similarity indicates the adequacy of the chosen multilayer skin model. The bright spots at the extrema of the sinusoidal boundary shape originate from strong back-reflection at the sections of layer boundaries normal to the incident probing beam direction.
Introduction of “speckle effects” into Monte Carlo simulation process with accordance to expression (8) causes additional noise observed in the obtained OCT images for all considered coherence length values (Fig. 4 ). These results are in qualitative agreement with the images simulated without account of the “speckle noise”, which on the one hand confirms the adequacy of the proposed algorithm and on the other hand demonstrates its ability for quantitative study of “speckle effects” in OCT image formation. In experimental setups “speckle effects” are usually reduced by spatial averaging of the signal and additional numerical processing.
At the next step we consider the implementation of expressions (16,17) into MC simulation which allows to account for polarization state of photons. The simulated A-scans (transversal position: x = 0 μm) for the multilayer skin model are presented in Fig. 5a . At small probing depth the co-polarization signal is extremely close to the non-polarized one while their discrepancy grows with the depth which is the evidence of depolarization of the probing radiation.
Figures 5b, 5c and 5d demonstrate the simulated OCT images obtained in non-, co- and cross- polarization modes correspondingly, the same as in Fig. 1 which represents similar experimental results. These images qualitatively agree with the experimental ones. Note that both experimental and simulated images obtained in cross-polarization mode demonstrate the same features: relatively low signal from the top boundary and bright epidermis layer. However, the simulated cross-polarization mode image exhibits lower absolute intensity compared to the experimental one which possibly originates from the assumption of soft scatterers accounted in expression (10). In reality the rate of depolarization appears to be higher compared to the value accounted in the simulations.
In this work, original theory for scattering of polarized radiation is developed and numerical simulations of OCT images are performed with implementation of this theory. Expressions for co- and cross-polarized radiation components contribution to OCT images of a biotissue sample are obtained. To simulate PS-OCT images of skin a sophisticated model mimicking physiological skin structure is embedded into MC code. Experimental OCT images of human skin are used as a reference. The parameters for skin model are chosen basing on the values available from literature.
We show that the developed improvement of MC code allows to simulate OCT images which are in good qualitative agreement with the experimental results. Simulation can be performed both with and without account for “speckle effects” which quantitatively demonstrates the role of these effects in OCT image formation for different objects. The abilities of the developed code are illustrated by good correlation with the OCT images of human skin in vivo.
The developed algorithm can be effectively utilized for simulation of polarization-sensitive OCT imaging of various scattering objects including biotissues and for evaluation of the role of various physical and instrumental factors such as speckles, temporal and/or spatial coherence properties of incident radiation, etc. Beside those, the phenomena of optical anisotropy of tissues such as birefringence and optical activity that are of a high interest in framework of modern biomedical optical diagnostics and can be introduced in the proposed algorithm with easier, e.g. similar to the approach presented in .
This project was supported in part by the Royal Society, UK, NATO (grant no. PST.CLG.979652), GETA graduate school and Tauno Tönning Foundation (Finland). E. Sergeeva, M. Kirillin and V. Kuzmin acknowledge Russian Foundation for Basic Research (grants 09-02- 97040, 10-02-00744, 10-02-00937) and grants of the President of Russian Federation МК-698.2009.2, МК-1127.2010.2. The work is partly supported by FTP “Scientific and scientific-pedagogical personnel of innovative Russia” (projects # 02.740.11.0839, 02.740.11.0566, 02.740.11.0437).
References and links
1. D. Huang, E. A. Swanson, C. P. Lin, J. S. Schuman, W. G. Stinson, W. Chang, M. R. Hee, T. Flotte, K. Gregory, C. A. Puliafito, and J. D. Fujimoto, “Optical coherence tomography,” Science 254(5035), 1178–1181 (1991). [CrossRef] [PubMed]
2. A. F. Fercher, W. Drexler, C. K. Hitzenberger, and T. Lasser, “Optical coherence tomography – principles and applications,” Rep. Prog. Phys. 66(2), 239–303 (2003). [CrossRef]
3. J. M. Schmitt, “Optical coherence tomography (OCT): A review,” IEEE J. Sel. Top. Quantum Electron. 5(4), 1205–1215 (1999). [CrossRef]
4. B. E. Bouma, and G. J. Tearney, Handbook of Optical Coherence Tomography, (Marcel Dekker, New York, 2002).
5. V. V. Tuchin, Handbook of Coherent Domain Optical Methods: Biomedical Diagnostics Environment and Material Science (Kluwer Academic, Boston, 2004).
6. M. J. Yadlowsky, J. M. Schmitt, and R. F. Bonner, “Multiple-scattering in optical coherence microscopy,” Appl. Opt. 43(25), 5699–5707 (1995). [CrossRef]
7. 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).
8. M. Yu. Kirillin, A. V. Priezzhev, and R. Myllylä, “Role of multiple scattering in formation of OCT skin images,” Quantum Electron. 38, 486–490 (2008). [CrossRef]
10. C. Lavigne, A. Roblin, V. Outters, S. Langlois, T. Girasole, and C. Roze, “Comparison of iterative and monte carlo methods for calculation of the Aureole about a point source in the earth’s atmosphere,” Appl. Opt. 38(30), 6237–6246 (1999). [CrossRef]
12. E. Berrocal, D. L. Sedarsky, M. E. 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]
13. 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–468 (2006). [CrossRef]
15. M. Yu. Kirillin, M. V. Shirmanova, M. A. Sirotkina, M. L. Bugrova, B. N. Khlebtsov, and E. V. Zagaynova, “Contrasting properties of gold nanoshells and titanium dioxide nanoparticles for OCT imaging of skin: Monte Carlo simulations and in vivo study,” J. Biomed. Opt. 14, 021017 (2009). [CrossRef] [PubMed]
16. B. Karamata, M. Laubscher, M. Leutenegger, S. Bourquin, T. Lasser, and P. Lambelet, “Multiple scattering in optical coherence tomography. I. Investigation and modeling,” J. Opt. Soc. Am. A 22(7), 1369–1379 (2005). [CrossRef]
17. B. Karamata, M. Leutenegger, M. Laubscher, S. Bourquin, T. Lasser, and P. Lambelet, “Multiple scattering in optical coherence tomography. II. Experimental and theoretical investigation of cross talk in wide-field optical coherence tomography,” J. Opt. Soc. Am. A 22(7), 1380–1388 (2005). [CrossRef]
18. M. Y. Kirillin, A. V. Priezzhev, and I. V. Meglinski, “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]
19. V. L. Kuzmin and I. V. Meglinski, “Multiple scattering and intensity fluctuations in optical coherence tomography of randomly inhomogeneous media,” J. Exp. Theor. Phys. 105(2), 285–291 (2007). [CrossRef]
20. R. V. Kuranov, V. V. Sapozhnikova, N. M. Shakhova, V. M. Gelikonov, E. V. Zagainova, and S. A. Petrova, “Combined application of optical methods to increase the information content of optical coherent tomography in diagnostics of neoplastic processes,” Quantum Electron. 32(11), 993–998 (2002). [CrossRef]
21. M. Yu. Kirillin, E. Alarousu, T. Fabritius, R. Myllylä, and A. V. Priezzhev, “Visualization of paper structure by optical coherence tomography: Monte Carlo simulations and experimental study,” J. Europ. Opt. Soc. Rap. Public. 2, 07031 (2007). [CrossRef]
22. I.M. Sobol’, The Monte Carlo Method (The University of Chicago Press, Chicago, 1974).
23. A. Ishimaru, Wave Propagation and Scattering in Random Media (Academic, New York, 1978).
25. D. Y. Churmakov, V. L. Kuz’min, and I. V. Meglinski, “Application of the vector Monte-Carlo method in polarisation optical coherence tomography,” Quantum Electron. 36(11), 1009–1015 (2006). [CrossRef]
26. J. W. Goodman, Statistical Optics (Wiley-Interscience, 1985).
27. C. Brosseau, Fundamentals of Polarized Light: a Statistical Optics Approach (New York: John Wiley & Sons, 1998).
28. C. F. Bohren, and D. R. Huffman, Absorption and scattering of light by small particles (New York: Wiley, 1983)
30. S. Bartel and A. H. Hielscher, “Monte Carlo simulations of the diffuse backscattering mueller matrix for highly scattering media,” Appl. Opt. 39(10), 1580–1588 (2000). [CrossRef]
31. M. J. Raković, G. W. Kattawar, M. B. Mehrubeoğlu, B. D. Cameron, L. V. Wang, S. Rastegar, and G. L. Coté, “Light backscattering polarization patterns from turbid media: theory and experiment,” Appl. Opt. 38(15), 3399–3408 (1999). [CrossRef]
32. D. A. Zimnyakov, Y. P. Sinichkin, P. V. Zakharov, and D. N. Agafonov, “Residual polarization of non-coherently backscattered linearly polarized light: the influence of the anisotropy parameter of the scattering medium,” Waves Random Media 11(4), 395–412 (2001). [CrossRef]
33. S. V. Gangnus, S. J. Matcher, and I. V. Meglinski, “Monte Carlo modeling of polarized light propagation in biological tissues,” Laser Phys. 14, 886–891 (2004).
34. J. M. Schmitt, A. H. Gandjbakhche, and R. F. Bonner, “Use of polarized light to discriminate short-path photons in a multiply scattering medium,” Appl. Opt. 31(30), 6535–6546 (1992). [CrossRef] [PubMed]
35. E. Akkermans, P. E. Wolf, R. Maynard, and G. Maret, “Theoretical-Study of the Coherent Backscattering of Light by Disordered Media,” J. Phys. France 49(1), 77–98 (1988). [CrossRef]
38. D. A. Zimnyakov and Y. P. Sinichkin, “A study of polarization decay as applied to improved imaging in scattering media,” J. Opt. A, Pure Appl. Opt. 2(3), 200–208 (2000). [CrossRef]
39. A. Dogariu, C. Kutsche, P. Likamwa, G. Boreman, and B. Moudgil, “Time-domain depolarization of waves retroreflected from dense colloidal media,” Opt. Lett. 22(9), 585–587 (1997). [CrossRef] [PubMed]
40. V.V. Tuchin, Tissue optics: light scattering methods and instruments for medical diagnosis (SPIE Press, Bellingham, 2000).
41. V. L. Kuzmin and I. V. Meglinski, “Helicity flip of backscattered circularly polarized light,” Proc. SPIE 7573, 75730Z (2010). [CrossRef]
42. P. S. Carney, E. Wolf, and G. S. Agarwal, “Statistical generalizations of the optical cross-section theorem with application to inverse scattering,” J. Opt. Soc. Am. A 14(12), 3366–3371 (1997). [CrossRef]
43. V. L. Kuzmin and E. V. Aksenova, “A generalized Milne solution for the correlation effects of multiple light scattering with polarization,” J. Exp. Theor. Phys. 96(5), 816–831 (2003). [CrossRef]
44. V. L. Kuz’min, V. P. Romanov, and L. A. Zubkov, “Propagation and scattering of light in fluctuating media,” Phys. Rep. 248(2-5), 71–368 (1994). [CrossRef]
45. P. K. Milsom, “A ray-optic, Monte Carlo, description of a Gaussian beam waist – applied to reverse saturable absorption,” Appl. Phys. B 70(4), 593–599 (2000). [CrossRef]
46. 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]