第1页 / 共10页
第2页 / 共10页
第3页 / 共10页
第4页 / 共10页
第5页 / 共10页
第6页 / 共10页
第7页 / 共10页
第8页 / 共10页
多层介质中偏振光传播的电场蒙特卡罗模拟.pdf
Electric field Monte Carlo simulation
of polarized light propagation in multi-
layered media
Chizhu Ding, Zuojun Tan, Shuhui Zhang, Siyu Chen
Chizhu Ding, Zuojun Tan, Shuhui Zhang, Siyu Chen, "Electric field Monte
Carlo simulation of polarized light propagation in multi-layered media," Proc.
SPIE 10461, AOPC 2017: Optical Spectroscopy and Imaging, 104610E (24
October 2017); doi: 10.1117/12.2282397
Event: Applied Optics and Photonics China (AOPC2017), 2017, Beijing, China
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 12/10/2017 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
PROCEEDINGS OF SPIESPIEDigitalLibrary.org/conference-proceedings-of-spie
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 12/10/2017 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
Electric field Monte Carlo simulation of polarized light propagation in multi-layered media Chizhu Dinga *, Zuojun Tana, Shuhui Zhanga, Siyu Chena aCollege of Science, Huazhong Agricultural University, Wuhan, Hubei, China ABSTRACT Electric field Monte Carlo (EMC) simulation is capable of modeling the polarization and coherence phenomena of light. Previous EMC program treat the turbid media as an infinite slab. An electric field Monte Carlo simulation of polarized light propagation in multi-layered media (EMCML) is presented in this paper. The complex electric field vectors are traced during the scattering and the reflection (or refraction) events. In order to improve the computational efficiency, our EMCML program is implemented in parallel in the GPU. The validity of EMCML is demonstrated by comparison between simulation results obtained by EMCML and previously reported programs. Keywords: Light propagation through tissue, photon migration, multiple scattering, backscattering, turbid media 1. INTRODUCTION Optical spectroscopy technique has been employed to noninvasively extract information from biological tissues. This information regarding biochemical composition and microarchitecture is useful for tissue characterization and diagnosis [1, 2]. Because biological tissue often consists of heterogeneous multi-layered structures, it would be extremely advantageous to be able to characterize the scattering and absorption properties of different layers from measurements of the emerging light. Several methods have been implemented to analyze multi-layered scattering media, such as time of flight [3, 4], diffusive wave [5, 6] and continuous wave [7]. Polarization techniques offer an alternative method and have the advantages of being relatively simple and inexpensive [8-12]. Multiple scattering depolarizes light and the initial polarization of incident light will be maintained only for a certain number of scattering events. Therefore, polarization gating can be used to selectively probing the properties of different depths within the layered media. Although the fundamental principle of polarization gating is straightforward, detailed quantification of the effect of optical properties on different polarization states has not been completely built. Analytically solving the radiative transfer equation of polarized light in turbid media is mathematically complex. Monte Carlo (MC) simulations of polarized light transport into scattering media have been developed by research groups. Most MC methods trace the Stokes vector of light in simulation and ignore the correlation effect of multiply scattered light [13, 14]. Min Xu proposed a MC algorithm (EMC) based on tracing the complex electric field instead of the Stokes vector [15]. It provides more detailed information about the wave properties than just polarization status and has been implemented in investigations of coherence, phase retardation, and interference phenomena [16-19]. Yaru Wang et al. optimized EMC by parallel computing using a graphics processing unit (GPU) under the compute unified device architecture (CUDA) and named it CUDAEMC [20]. However, both EMC and CUDAEMC model the turbid media as a single homogeneous slab and are not suitable for multi-layered media. We report a GPU accelerated electric field MC simulation of polarized light propagation in multi-layered media, called EMCML. It is capable of simulating the coherence and polarization phenomena of light within multi-layered media. The detailed implementation of EMCML is described in Section 2. EMCML computation results are presented and analyzed in Section 3. The conclusion is given in Section 4. *dingchizhu@mail.hzau.edu.cn AOPC 2017: Optical Spectroscopy and Imaging, edited by Jin Yu, Zhe Wang, Wei Hang, Bing Zhao, Xiandeng Hou, Mengxia Xie, Tsutomu Shimura, Proc. of SPIE Vol. 10461, 104610E © 2017 SPIE · CCC code: 0277-786X/17/$18 · doi: 10.1117/12.2282397Proc. of SPIE Vol. 10461 104610E-1
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 12/10/2017 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
2. METHODS 2.1 Structure of the simulation MC simulation code of unpolarized light transport in multi-layered tissues (MCML) has been developed and brought to public domain by Wang and Jacques [21]. We incorporated the polarization effect into the MCML code by tracing the parallel and perpendicular components of the complex electric field with respect to the scattering plane or the incidence plane. The main steps of photon tracing are briefly summarized below. 1. Launch a photon packet with initial photon direction, position, weight and parallel and perpendicular components of the complex electric field. 2. Generate a random step size based on a negative exponential distribution with the mean length ls = 1/μs, where μs is the scattering coefficient. 3. Check whether the photon will hit the boundary or interface. If not, move the photon, attenuate the photon weight, scatter the photon (will be explained in details later), and then skip to Step 5; otherwise, proceed to Step 4. 4. Move the photon to the boundary or interface. Determine whether the photon is transmitted or reflected (will be explained in details later), and then transmit or reflect the photon. 5. Perform photon survival check. If the photon exits the media or the number of scattering events is larger than a predefined number, save the temporary simulation results and end the simulation for the photon; if the weight of the photon is lower than a predefined threshold value, perform roulette to update the weight of the photon; otherwise, go back to Step 2. 6. Repeat from Step 1 until all the photon simulations have been executed. 2.2 Scattering events The EMCML program updates the complex electric field vectors during scattering events in the same way as that in the original EMC code for a homogeneous semi-infinite medium. A detailed introduction to the EMC can be found in the literature [15]. To improve understanding of scattering events, key points are described below. As depicted in Figure 1(a), s and s’ are the unit vectors in the directions of the photon propagation before and after a scattering event, respectively. The current scattering plane is the plane spanned by s and s’. The local orthonormal coordinate system (m, n, s) is built where m and n are the unit vectors in the directions of the parallel and perpendicular components E1 and E2 of the electric field with respect to the scattering plane of the previous scattering event. For scattering with the scattering angle θ and the azimuthal angle φ, the local coordinate system (m, n, s) is rotated to (m’, n’, s’) as coscoscossinsinsincos0sincossinsincosθϕθϕθϕϕθϕθϕθ′−⎛⎞⎛⎞⎛⎞⎜⎟⎜⎟⎜⎟′=−⎜⎟⎜⎟⎜⎟⎜⎟⎜⎟⎜⎟′⎝⎠⎝⎠⎝⎠mmnnss (1) If the medium is composed of Mie scatterers as in most biological tissues [22], then 2×2 the amplitude scattering matrix S(θ) can be derived from Mie theory. The parallel and perpendicular components of the electric field are scattered according to Sj(θ) where j = 2, 1, respectively. The electric field E = E1m + E2n is scattered to E’ = E1’m’ + E2’n’ by []1/212212112cossin(,)sincosESSEFESSEϕϕθϕϕϕ−′⎛⎞⎛⎞⎛⎞=⎜⎟⎜⎟⎜⎟′−⎝⎠⎝⎠⎝⎠ (2) with a factor Proc. of SPIE Vol. 10461 104610E-2
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 12/10/2017 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
INeli=m'(a)(b) ()()()()222222222221121222*2112(,)cossinsincos2cossinReFSSESSESSEEθϕϕϕϕϕϕϕ=+++⎡⎤+−⎣⎦ (3) to normalize the light intensity. Figure 1 Local coordinate system transformation: (a) for a scattering event; (b) for a reflection or refraction event at the boundary or interface The probability density function of θ can be calculated as ()()()()2221220,dscaSSppxQθθθθϕϕπ+==∫ (4) where x is the size parameter and x = 2πna/λ; n is the refractive index of the host medium, a is the particle radius and λ is the wavelength of light in the host medium. The source code developed by Bohren and Huffman [23] is used to calculate the amplitude scattering matrix element S1(θ) and S2(θ). The rejection technique [7] is used to sample the scattering angle θ and the azimuthal angle φ. 2.3 Reflection or refraction events When a photon hits a boundary of the tissue or an interface between layers in the tissue with a mismatched refractive index, it is either reflected or transmitted. Two steps of local coordinate system transformations are performed as depicted in Figure 1(b). 1. The local coordinate system (m, n, s) is rotated about s until m enters into the plane of incidence, and n is perpendicular to the plane of incidence. The two components of the electric field are updated to guarantee E1m + E2n = E1’m’ + E2’n’. 2. The internal reflectance R is calculated according to Fresnel’s law [24]. A random number ξ is generated and compared with R to determine whether the photon will be transmitted or reflected. 3. The unit vector st (or sr) in the direction of the transmitted (or reflected) photon is determined by Snell’s law (or laws of reflection); the perpendicular vector nt (or nr) remains the same as n’; and the parallel vector mt = nt × st (or mr = nr × sr). 2.4 GPU implementation MC methods are statistics-based computational algorithms that need to trace a large number of independent photons to achieve accurate results. Therefore, conventional MC simulations are time-consuming. Parallel computing with GPU has Proc. of SPIE Vol. 10461 104610E-3
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 12/10/2017 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
been utilized to accelerate MC simulations by several research groups [20, 25-27]. We paralyzed the computation based on CUDA published by NVIDIA Corporation. When the main program starts in the CPU, tables of S1 and S2 (the elements of the amplitude scattering matrix) for each layer are generated by Mie theory. They will be copied to the GPU memory and be used to calculate the scattering of photons later in the GPU. Then, the GPU kernel is started and each thread traces one photon packet independently. The number of the photon packets being traced simultaneously depends on the number of threads used by the GPU. Because there is generally a significant number of GPU threads, the computational time is greatly less than that required by the conventional MC simulation in which the photons are traced one by one. After all the photons have been completed, the simulation results are copied back to the CPU and saved. 3. RESULTS AND DISCUSSION In order to validate our EMCML program, we compare the simulation results by our program with the results by the reported programs. Two kinds of tissue models are used. One is an infinite slab model, the other is a three-layer model as shown in Figure 2. Parameters of the tissue model are listed in Table 1. The tissues are assumed to be made up of uniformly distributed small particles and the optical parameters of the tissues are obtained by Mie theory. The optical parameters of the first, the second and the third layers are of similar values with those in the epidermis, the dermis and the blood plexus layers of skin model at 442 nm [28], respectively. The incident light propagates in the z direction, normally incident on the surface of the tissue model. The wavelength of the incident light is 442 nm and a total of 108 photons are launched in each simulation. Figure 2 Illustration of the three-layer model used in the EMCML simulation Table 1 Parameters of the tissue model for 442 nm laser Tissue Type Layer 1 Layer 2 Layer 3 nmed 1.4 1.4 1.37 npar 1.59 1.59 1.5 radius (μm) 0.156 0.156 4.5 ρ (μm-3) 2.332 2.865 0.00035 μa (cm-1) 56 6.7 541 μs (cm-1) 569.9 700.2 521.7 g 0.8 0.8 0.95 d (cm) 1-layer model 0.03 -- -- 3-layer model 0.005 0.015 0.01 3.1 Polarized light propagation in an infinite slab model The first validation is done by simulating polarized light propagation in the single-layer model described in Table 1. Three kinds of polarization states of the incident light are computed. They are the horizontal polarized light (polarized in the x direction), the 45-degree polarized light and the right circular polarized light. The difference between our EMCML program and Min Xu’s EMC [15] results is always less than 0.2%. Proc. of SPIE Vol. 10461 104610E-4
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 12/10/2017 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
0.01IE 4- MCML- EMCML- EMCML with HG0.0000.0050.0100.0150.020source -detector seperation [cm]Lio:81000.010.000.010.02depth [cm]0.03 We also compare the run time by EMC and EMCML programs and list the run time in Table 2. The EMC program was run using an Intel Xeon E5-2620 v3 CPU, while the EMCML program used a NVIDIA Quadro K2200 GPU. For the above simulations, a speedup of over 45x is achieved by the parallel computing. More simulations with different parameters demonstrate that the parallel computing speed would be affected by the number of program branches in the GPU, which is related to the mean scattering length and the thickness of the tissue model. Table 2 Comparison of the run time by EMC and EMCML Polarization of incident light EMC EMCML Speedup Horizontal polarized 15441 340.8 45.3 45-degree polarized 15403 340.9 45.2 Right circular polarized 15391 340.7 45.2 3.2 Unpolarized light propagation in a multi-layered model The second validation is carried out by simulating unpolarized light propagation in the three-layer model described in Table 1. The incident light is modeled as randomly polarized incident light in our EMCML program. And the results are compared with those obtained by Wang and Jacques’s MCML program [21] without taking into account the polarization effects. The spatial distribution of diffuse reflectance and energy deposition are plotted in Figure 3(a) and (b). The energy deposition curves in Figure 3(b) are divided into three sections because the tissue model comprises three layers and each layer has a different absorption coefficient. There are differences between the curves obtained by MCML and EMCML programs. Some numerical results are summarized in Table 3. A maximum relative error of over 5% is found between the results. It might be because a different sampling strategy (the rejection technique [7]) is used to sample the scattering angle θ and the azimuthal angle φ in our EMCML program. As can be seen in Figure 3 and Table 3, if we replace the rejection method in EMCML by the conventional sampling strategy according to the Henyey-Greenstein (HG) function, good agreement is observed between the simulation outputs obtained by MCML and EMCML with HG function. (a) (b) Figure 3 Comparison of the run results by MCML, EMCML and EMCML with HG function. (a) Diffuse reflectance; (b) Absorption. Table 3 Comparison of the run results by MCML, EMCML and EMCML with HG function MCML EMCML EMCML with HG Result Error (%) Result Error (%) Specular reflectance 0.0277778 0.0277778 0 0.0277778 0 Diffuse reflectance 0.126887 0.125579 -1.04 0.126953 0.052 Absorption 0.845158 0.846501 0.16 0.845133 -0.003 by layer 1 0.5151 0.521104 1.15 0.515639 0.105 by layer 2 0.1185 0.118701 0.17 0.118275 -0.190 by layer 3 0.2116 0.206696 -2.37 0.21122 -0.180 Transmittance 0.00017735 0.00018677 5.05 0.00017815 0.453 Proc. of SPIE Vol. 10461 104610E-5
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 12/10/2017 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
1 +1-7-6.5-6-5.5-5-4.5-4-3.5DR00.20.40.60.81 3.3 Polarized light propagation in a multi-layered model One advantage of our EMCML program is the ability to deal with polarized light propagation in a multi-layered model. We simulate light of different initial polarization states propagating through the three-layer tissue model in Table 1. The average backscattering light intensity Ix + Iy, the x component Ix, the y component Iy and the depolarization ratio are displayed as two-dimensional images in Figure 4Error! Reference source not found.. The depolarization ratio DR of backscattered light intensity is defined as: xyxyIIDRII−=+ (5) Images produced by the horizontal polarized light are presented. Each computed image is made up of a 100 × 100 detection grid, with the laser incident in the center. The side length is 20ls1. Figure 4 Backscattering light intensity and depolarization ratio produced by horizontal polarized incident light. The first three images display the average backscattering light intensity Ix + Iy, the x component Ix and the y component Iy, respectively. The fourth image displays the distribution of depolarization ratio DR. In EMCML program, we record the layer the photons penetrate and distinguish the x and y components of the backscattering signals generated by photons penetrating into different layers respectively. The cross sections of the distribution of light intensity along the x axis are plotted in Figure 5. For both the first and the second layers, the co-polarized components Ix are higher than the cross-polarized components Iy. That is because the co-polarized signal Ix is comprised of both singly scattered photons primarily from shallow depths and multiply scattered photons from deeper depths. And the cross-polarized signal Iy is mainly comprised of multiply scattered photons [10]. The light intensity generated by photons penetrating into the third layer is so weak that it is ignorable. The differential polarization signals Ix – Iy along the x axis generated by photons from the first and the second layers are depicted in Figure 6. In the source-detector separation less than 5ls1, the differential polarization signal is predominantly determined by the photons from the first layer. Therefore, polarization gating technique helps selectively probe the properties of superficial tissue. Proc. of SPIE Vol. 10461 104610E-6
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 12/10/2017 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
- Layer 1, Ix-- Layer 1, /y1E -3-Layer 2,Ix-- Layer 2, /y- Layer 3, Ix-- Layer 3, /yN 1E-4.lE-8o5t' [to]Io1Er3=1Er4-1E-51E-6-1E-7-1E 8- Layerl- Layer2510 Figure 5 The x and y components of backscattering light intensity along the x axis generated by photons from different layers. Figure 6 The differential polarization signals Ix – Iy along the x aixs generated by photons from layer 1 and layer 2. 4. CONCLUSION In this paper, we present a GPU accelerated electric field Monte Carlo simulation of polarized light propagation in multi-layered media (EMCML). Because biological tissue often consists of heterogeneous multi-layered structures, the EMCML program would be more suitable for analyze the polarization and coherence phenomena in biological optics than the previous EMC program. ACKNOWLEDGEMENTS The research was financially supported by the National Natural Science Foundation of China (Grant No. 61505063) and the Fundamental Research Funds for the Central Universities (No. 2662015BQ046). Proc. of SPIE Vol. 10461 104610E-7
相关推荐
2023年江西萍乡中考道德与法治真题及答案.doc
2012年重庆南川中考生物真题及答案.doc
2013年江西师范大学地理学综合及文艺理论基础考研真题.doc
2020年四川甘孜小升初语文真题及答案I卷.doc
2020年注册岩土工程师专业基础考试真题及答案.doc
2023-2024学年福建省厦门市九年级上学期数学月考试题及答案.doc
2021-2022学年辽宁省沈阳市大东区九年级上学期语文期末试题及答案.doc
2022-2023学年北京东城区初三第一学期物理期末试卷及答案.doc
2018上半年江西教师资格初中地理学科知识与教学能力真题及答案.doc
2012年河北国家公务员申论考试真题及答案-省级.doc
2020-2021学年江苏省扬州市江都区邵樊片九年级上学期数学第一次质量检测试题及答案.doc
2022下半年黑龙江教师资格证中学综合素质真题及答案.doc
