<\/p>\n\n\n\n
We store the information we want to keep (in our case, the reflectance) in \u201cbins\u201d that we add to at each iteration stage.<\/p>\n\n\n\n Here\u2019s the general flow for Monte Carlo simulation:<\/p>\n\n\n\n To translate this: We launch <\/em><\/strong>a photon with a weight (w<\/em> = 1.0), we move the photon (hop<\/strong><\/em>) and check if we have changed boundaries between tissues. At the drop<\/em><\/strong> stage we interact with the tissue (absorption, etc). Next, we redirect the photon based on scattering (spin<\/em><\/strong>), and finally we run a simple roulette to see if the photon should terminate<\/em><\/strong>. If we don\u2019t terminate, we keep moving the photon. If it does <\/em>terminate, we launch the next photon. This process continues until we run out of photons (N photons<\/em><\/strong>).<\/p>\n\n\n\n <\/p>\n\n\n\n We will now cover each of these steps more deeply:<\/p>\n\n\n\n To launch we need three things: The first two should be fairly explanatory. We start at (0,0,0) and initially move forwards \u2013 we will use +z<\/em> as our forward axis into the tissue. <\/p>\n\n\n\n <\/p>\n\n\n\n We will use an \u201cisotropic point source\u201d to launch our photon. This simply means we have no preferred direction. Our starting position shall be:<\/p>\n\n\n\n x = 0, y = 0, z = 0<\/p>\n\n\n\n And our launching trajectory will be:<\/p>\n\n\n\n RND = 0.0 <= n<\/em> < 1.0 <\/p>\n\n\n\n ux = sin\u03b8 cos\u03d5 <\/p>\n\n\n\n The weight, w<\/em>, is what we shall use to track how much of the photon\u2019s energy is remaining. In our simulation we shall define a zone around the origin through which we will collect any photons that head back <\/em>into it. By summing the weight of these photons, at all the wavelengths of light we are sampling, we will end up with our reflectance values.<\/p>\n\n\n\n <\/em><\/p>\n\n\n\n Hop consists of two parts: A Boundary Check<\/em> and Moving<\/em>. In our two-layer model, we have an upper epidermis, a lower dermis layer and the external medium \u2013 air. <\/p>\n\n\n\n <\/p>\n\n\n\n For our purposes, we count anything above the epidermis as \u2018escaped\u2019, and if it is close enough to the beam source (the launch<\/strong><\/em> position), we capture its weight.<\/p>\n\n\n\n <\/p>\n\n\n\n We\u2019ll start by calculating our new position, this is calculated by:<\/p>\n\n\n\n \u00b5t<\/sub> = \u00b5a<\/sub> + \u00b5s<\/sub> <\/p>\n\n\n\n s = -ln(RND)\/\u00b5t<\/sub> <\/p>\n\n\n\n Now we check whether it has escaped the tissue by checking if the uz property is < 0.0. However<\/span>, <\/em>the photon may partially reflect back into the tissue at the boundary, due to total internal reflection. <\/p>\n\n\n\n <\/p>\n\n\n\n How do we calculate this?<\/p>\n\n\n\n First, let\u2019s start by moving the photon back to the point where it left the tissue, and move it just back to the surface. We calculate the partial step, s1<\/sub>:<\/p>\n\n\n\n <\/p>\n\n\n\n s1<\/sub> = abs(z\/uz)<\/p>\n\n\n\n <\/p>\n\n\n\n Then, retract the full step:<\/p>\n\n\n\n x = x \u2013 s ux and then add in the partial step, back to the boundary:<\/p>\n\n\n\n x = x + s1<\/sub> ux <\/p>\n\n\n\n Now our photon is at the surface, we use the Fresnel equation to determine how much of the energy was reflected back into the tissue:<\/p>\n\n\n\n We can now calculate the weight that has escaped the tissue with: <\/p>\n\n\n\n Re<\/sub> = 1 \u2013 Ri<\/sub><\/p>\n\n\n\n <\/p>\n\n\n\n This can then be stored in our \u2018bin\u2019<\/p>\n\n\n\n r = \u221a<\/strong>(x2<\/sup> + y2<\/sup>) <\/p>\n\n\n\n where dr is our radial bin size, (radial size of our photons\/number of bins). I have glossed over these parameters, as its a bit deep for this tutorial. Simply put, r is the radial position of the photon, dr gives us a \u2018depth\u2019. Although we simulate the Monte Carlo simulation in 3D Cartesian coordinates, we actually only care about getting data for radial position and depth, as the photon beam is cylindrical.<\/p>\n\n\n\n Finally, we bounce the photon back into the tissue, with it\u2019s escapes weighting removed:<\/p>\n\n\n\n w = Ri<\/sub> * w <\/p>\n\n\n\n Our drop stage is the easiest. As we have calculated the reflectance of the tissue above, we simply need to decrement the weight of our photon by the amount absorbed by the surface albedo.<\/p>\n\n\n\n Which albedo value we use depends on the layer in the skin we are in:<\/p>\n\n\n\n if z <= epidermis_thickness <\/p>\n\n\n\n Now to reduce the weight:<\/p>\n\n\n\n absorb = w * (1 \u2013 albedo) Next up: scatter the photon! We need two angles of scatter, \u03b8 (deflection) and \u03d5 (azimuthal). The source paper goes into great detail as to why<\/em> we use the following formula (the HG function) to calculate \u03b8, if you want more information, but for the purposes of this tutorial, we\u2019ll just describe it:<\/p>\n\n\n\n The azimuthal angle is a bit less verbose:<\/p>\n\n\n\n \u03d5 = 2\u03c0 RND<\/p>\n\n\n\n Now we can calculate the new angles, we need to update our photon trajectories:<\/p>\n\n\n\n if (1 \u2013 abs(uz) <= 1 \u2013 cos0) : # close to perpendicular. ux = uxx Where SIGN(x) is x >= 0 ? 1 : 0<\/p>\n\n\n\n <\/p>\n\n\n\n The photon will continue to propagate naturally and will never terminate naturally (w<\/em> will simply get smaller). We therefore need a way to terminate the photon. <\/p>\n\n\n\n This is done by simply creating a random number and comparing it to the remaining energy of the photon, if the photon weight is below a threshold:<\/p>\n\n\n\n if w < THRESHOLD However, this will not conserve energy as we are removing [the leftover] energy from the simulation (and energy cannot be destroyed). We therefore need to conserve the energy.<\/p>\n\n\n\n We do this by adding energy into photons that are not terminated by the random chance, i.e., if we expect that 9\/10 photons will terminate with the random check, we multiply the energy by 10-fold on the one photon that did not terminate and allow it to continue. This maintains energy conservation laws. Therefore, the actual implementation is as follows:<\/p>\n\n\n\n if w < THRESHOLD The final step in the simulation, after letting it run through as many iterations as we want, is to get the data back out.<\/p>\n\n\n\n To do this, we simply iterate over the \u2018bins\u2019 we added weights into earlier, summing up their weights, and diving it down by the number of photons we used in our simulation:<\/p>\n\n\n\n for i = 0; i < length(reflectionBin); i++ Where Rt<\/sub> is total reflection<\/p>\n\n\n\n <\/p>\n\n\n\n Now we do this for each wavelength we want to sample, and we will get a nice set of data describing the reflection at each of those points. We\u2019ll implement this all in Python in the next part.<\/p>\n","protected":false},"excerpt":{"rendered":" Introduction As discussed in previous chapters, we need to simulate light photon beams in order to determine the true color of our surface. The \u201cMonte Carlo\u201d simulation provides a method to do this. The name \u201cMonte Carlo\u201d derives from the randomness used in the method; we use random values at several stages and also use a Russian Roulette approach to \u2018killing\u2019 the photon beam. The general approach is to create a beam of light at a position (i.e. the origin, (0,0,0)), give it a direction and wavelength. We then iterate for n times, each time moving the photon beam further into the material, scattering, absorbing and reflecting as it goes. We run the simulation for each wavelength of light we care about \u2013 i.e. 380nm-780nm, at 10nm intervals. We store the information we want to keep (in our case, the reflectance) in \u201cbins\u201d that we add to at each iteration stage. Here\u2019s the general flow for Monte Carlo simulation: To translate this: We launch a photon with a weight (w = 1.0), we move the photon (hop) and check if we have changed boundaries between tissues. At the drop stage we interact with the tissue (absorption, etc). Next, we redirect the photon based on scattering (spin), and finally we run a simple roulette to see if the photon should terminate. If we don\u2019t terminate, we keep moving the photon. If it does terminate, we launch the next photon. This process continues until we run out of photons (N photons). We will now cover each of these steps more deeply: Launch To launch we need three things:1. A starting location2. A starting direction to hop3. A \u201cweight\u201d The first two should be fairly explanatory. We start at (0,0,0) and initially move forwards \u2013 we will use +z as our forward axis into the tissue. We will use an \u201cisotropic point source\u201d to launch our photon. This simply means we have no preferred direction. Our starting position shall be: x = 0, y = 0, z = 0 And our launching trajectory will be: RND = 0.0 <= n < 1.0cos\u03b8 = 2 RND \u2013 1 sin\u03b8 = \u221a(1-cos\u03b8)\u03d5 = 2 \u03c0 RNDif ( \u03d5 < \u03c0 )\u2003sin\u03d5 = \u221a(1-cos2\u03d5)else\u2003sin\u03d5 = \u2013\u221a(1-cos2\u03d5) ux = sin\u03b8 cos\u03d5uy = sin\u03b8 sin\u03d5uz = abs(cos\u03b8) The weight, w, is what we shall use to track how much of the photon\u2019s energy is remaining. In our simulation we shall define a zone around the origin through which we will collect any photons that head back into it. By summing the weight of these photons, at all the wavelengths of light we are sampling, we will end up with our reflectance values. Hop Hop consists of two parts: A Boundary Check and Moving. In our two-layer model, we have an upper epidermis, a lower dermis layer and the external medium \u2013 air. For our purposes, we count anything above the epidermis as \u2018escaped\u2019, and if it is close enough to the beam source (the launch position), we capture its weight. We\u2019ll start by calculating our new position, this is calculated by: \u00b5t = \u00b5a + \u00b5swhere \u00b5a is the absorption coefficient, and \u00b5s is our scattering coefficient s = -ln(RND)\/\u00b5tx = x + s uxy = y + s uyz = z + s uz Now we check whether it has escaped the tissue by checking if the uz property is < 0.0. However, the photon may partially reflect back into the tissue at the boundary, due to total internal reflection. How do we calculate this? First, let\u2019s start by moving the photon back to the point where it left the tissue, and move it just back to the surface. We calculate the partial step, s1: s1 = abs(z\/uz) Then, retract the full step: x = x \u2013 s uxy = y \u2013 s uyz = z \u2013 s uz and then add in the partial step, back to the boundary: x = x + s1 uxy = y + s1 uyz = z + s1 uz Now our photon is at the surface, we use the Fresnel equation to determine how much of the energy was reflected back into the tissue: We can now calculate the weight that has escaped the tissue with: Re = 1 \u2013 Ri This can then be stored in our \u2018bin\u2019 r = \u221a(x2 + y2)ir = r\/drreflectionBin[ir] = Re where dr is our radial bin size, (radial size of our photons\/number of bins). I have glossed over these parameters, as its a bit deep for this tutorial. Simply put, r is the radial position of the photon, dr gives us a \u2018depth\u2019. Although we simulate the Monte Carlo simulation in 3D Cartesian coordinates, we actually only care about getting data for radial position and depth, as the photon beam is cylindrical. Finally, we bounce the photon back into the tissue, with it\u2019s escapes weighting removed: w = Ri * wuz = -uzx = (s-s1) * uxy = (s-s1) * uyz = (s-s1) * uz Drop Our drop stage is the easiest. As we have calculated the reflectance of the tissue above, we simply need to decrement the weight of our photon by the amount absorbed by the surface albedo. Which albedo value we use depends on the layer in the skin we are in: if z <= epidermis_thickness\u2003albedo = \u03c3epidermiss\/(\u03c3epidermiss + \u03c3epidermisa)else\u2003albedo = \u03c3dermiss\/(\u03c3dermiss + \u03c3dermisa) Now to reduce the weight: absorb = w * (1 \u2013 albedo)w = w \u2013 absorb Spin Next up: scatter the photon! We need two angles of scatter, \u03b8 (deflection) and \u03d5 (azimuthal). The source paper goes into great detail as to why we use the following formula (the HG function) to calculate \u03b8, if you want more information, but for the purposes of this tutorial, we\u2019ll just describe it: The azimuthal angle is a bit less verbose: \u03d5 = 2\u03c0 RND Now we can calculate the new angles, we need to update our photon trajectories: if (1 \u2013 abs(uz) <= 1 \u2013 cos0)<\/p>\n","protected":false},"author":1,"featured_media":254,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[1],"tags":[8,7],"class_list":["post-146","post","type-post","status-publish","format-standard","has-post-thumbnail","hentry","category-uncategorized","tag-lut","tag-skin"],"_links":{"self":[{"href":"https:\/\/half4.xyz\/index.php\/wp-json\/wp\/v2\/posts\/146"}],"collection":[{"href":"https:\/\/half4.xyz\/index.php\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/half4.xyz\/index.php\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/half4.xyz\/index.php\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/half4.xyz\/index.php\/wp-json\/wp\/v2\/comments?post=146"}],"version-history":[{"count":26,"href":"https:\/\/half4.xyz\/index.php\/wp-json\/wp\/v2\/posts\/146\/revisions"}],"predecessor-version":[{"id":1284,"href":"https:\/\/half4.xyz\/index.php\/wp-json\/wp\/v2\/posts\/146\/revisions\/1284"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/half4.xyz\/index.php\/wp-json\/wp\/v2\/media\/254"}],"wp:attachment":[{"href":"https:\/\/half4.xyz\/index.php\/wp-json\/wp\/v2\/media?parent=146"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/half4.xyz\/index.php\/wp-json\/wp\/v2\/categories?post=146"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/half4.xyz\/index.php\/wp-json\/wp\/v2\/tags?post=146"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}![\"\"](\"https:\/\/half4.xyz\/wp-content\/uploads\/2020\/07\/image-31.png\")
![\"\"](\"https:\/\/half4.xyz\/wp-content\/uploads\/2020\/07\/image-32.png\")
Launch<\/h2>\n\n\n\n
1. A starting location
2. A starting direction to hop
3. A \u201cweight\u201d<\/p>\n\n\n\n
cos\u03b8 = 2 RND \u2013 1
sin\u03b8 = \u221a<\/strong>(1-cos\u03b8)
\u03d5 = 2 \u03c0 RND
if ( \u03d5 < \u03c0 )
\u2003sin\u03d5 = \u221a<\/strong>(1-cos2<\/sup>\u03d5)
else
\u2003sin\u03d5 = \u2013\u221a<\/strong>(1-cos2<\/sup>\u03d5)<\/p>\n\n\n\n
uy = sin\u03b8 sin\u03d5
uz = abs(cos\u03b8)<\/p>\n\n\n\nHop<\/h2>\n\n\n\n
where \u00b5a<\/sub> is the absorption coefficient, and \u00b5s<\/sub> is our scattering coefficient<\/p>\n\n\n\n
x = x + s ux
y = y + s uy
z = z + s uz<\/p>\n\n\n\n
y = y \u2013 s uy
z = z \u2013 s uz<\/p>\n\n\n\n
y = y + s1<\/sub> uy
z = z + s1<\/sub> uz<\/p>\n\n\n\n![\"\"](\"https:\/\/half4.xyz\/wp-content\/uploads\/2020\/07\/image-34.png\")
ir = r\/dr
reflectionBin[ir] = Re<\/sub><\/p>\n\n\n\n
uz = -uz
x = (s-s1<\/sub>) * ux
y = (s-s1<\/sub>) * uy
z = (s-s1<\/sub>) * uz<\/p>\n\n\n\nDrop<\/h2>\n\n\n\n
\u2003albedo = \u03c3epidermis<\/sup>s<\/sub>\/(\u03c3epidermis<\/sup>s<\/sub> + \u03c3epidermis<\/sup>a<\/sub>)
else
\u2003albedo = \u03c3dermis<\/sup>s<\/sub>\/(\u03c3dermis<\/sup>s<\/sub> + \u03c3dermis<\/sup>a<\/sub>)<\/p>\n\n\n\n
w = w \u2013 absorb<\/p>\n\n\n\nSpin<\/h2>\n\n\n\n
![\"\"](\"https:\/\/half4.xyz\/wp-content\/uploads\/2020\/07\/image-35.png\")
\u2003uxx = sin\u03b8 * cos\u03d5
\u2003uyy = sin\u03b8 * sin\u03d5
\u2003uzz = cos\u03b8 * SIGN(uz)
else:
\u2003temp = \u221a<\/strong>(1.0 \u2013 uz * uz)
\u2003uxx = sin\u03b8 * (ux * uz * cos\u03d5 \u2013 uy * sin\u03d5) \/ temp + ux * cos\u03b8
\u2003uyy = sin\u03b8 * (uy * uz * cos\u03d5 + ux * sin\u03d5) \/ temp + uy * cos\u03b8
\u2003uzz = -sin\u03b8 * cos\u03d5 * temp + uz * cos\u03b8<\/p>\n\n\n\n
uy = uyy
uz = uzz<\/p>\n\n\n\nTerminate<\/h2>\n\n\n\n
\u2003if RND > CHANCE
\u2003\u2003terminate the photon<\/p>\n\n\n\n
\u2003if RND > CHANCE
\u2003\u2003terminate the photon
\u2003else
\u2003\u2003w = w\/CHANCE<\/p>\n\n\n\nExtracting The Reflectance<\/h2>\n\n\n\n
\u2003Rt<\/sub> = Rt<\/sub> + reflectionBin[i] \/Nphotons<\/p>\n\n\n\n