import random\nimport math\nclass SkinAbsorption:\n ....<\/code><\/pre>\n\n\n\nIn our Python script, we\u2019ll add a new function to our class. We\u2019ll call this MonteCarlo. As its inputs, we need these parameters from our spectral work: \u03c3epidermis<\/sup>a<\/sub>, \u03c3epidermis<\/sup>s<\/sub>, \u03c3dermis<\/sup>a<\/sub>, \u03c3dermis<\/sup>s<\/sub> and \u03bb<\/p>\n\n\n\ndef MonteCarlo (self, epi_mua, epi_mus, derm_mua, derm_mus, nm):<\/code><\/pre>\n\n\n\nWe\u2019ll need two more utility functions. One is the Fresnel function, the second is a quick sign(x) function. Add these to the class, too:<\/p>\n\n\n\n
def RFresnel(self, n1, n2, cosT1):\n r = 0.0\n cosT2 = 0.0\n COSZERO = 1.0 - 1.0e-12\n COS90D = 1 * pow(10,-6)\n if n1 == n2: #matched boundary\n r = 0.0\n cosT2 = cosT1\n elif cosT1 > COSZERO: # normal incident\n cosT2 = ca1\n r = (n2-n1)\/(n2+n1)\n r *= r\n elif cosT1 < COS90D: # very slant\n cosT2 = 0.0\n r = 1.0\n else: #general\n sinT1 = math.sqrt(1 - cosT1*cosT1)\n sinT2 = n1 * sinT1\/n2\n \n if sinT2 >= 1.0:\n r = 1.0\n cosT2 = 0.0\n else:\n cosT2 = math.sqrt(1 - sinT2 * sinT2)\n cosAP = cosT1*cosT2 - sinT1*sinT2\n cosAM = cosT1*cosT2 + sinT1*sinT2\n sinAP = sinT1*cosT2 + cosT1*sinT2\n sinAM = sinT1*cosT2 - cosT1*sinT2\n r = 0.5 * sinAM * sinAM*(cosAM*cosAM+cosAP*cosAP)\/(sinAP*sinAP*cosAM*cosAM)\n return r<\/code><\/pre>\n\n\n\nThe next step is to add the following values to our Monte Carlo function. These are values we need for our light transport:<\/p>\n\n\n\n
# These are our Monte Carlo Light Transport Variables that don't change\nNbins = 1000\nNbinsp1 = 1001\nPI = 3.1415926\nLIGHTSPEED = exp(2.997925,10)\nALIVE = 1\nDEAD = 0\nTHRESHOLD = 0.01\nCHANCE = 0.1\nCOS90D = exp(1,-6)\nONE_MINUS_COSZERO = exp(1,-12)\nCOSZERO = 1.0 - 1.0e-12 # cosine of about 1e-6 rad\ng = 0.9\nnt = 1.4 # Index of refraction\nepidermis_thickness = 0.25<\/code><\/pre>\n\n\n\nLets now create all the variables we need to change, such as position, trajectory, bins, weight, etc.:<\/p>\n\n\n\n
x = 0.0\ny = 0.0\nz = 0.0 # photon position\n\nux = 0.0\nuy = 0.0\nuz = 0.0 # photon trajectory as cosines\n\nuxx = 0.0\nuyy = 0.0\nuzz = 0.0 # temporary values used during SPIN\n\ns = 0.0 # step sizes. s = -log(RND)\/mus [cm] \ncostheta = 0.0 # cos(theta) \nsintheta = 0.0 # sin(theta) \ncospsi = 0.0 # cos(psi) \nsinpsi = 0.0 # sin(psi) \npsi = 0.0 # azimuthal angle \ni_photon = 0.0 # current photon\nW = 0.0 # photon weight \nabsorb = 0.0 # weighted deposited in a step due to absorption \nphoton_status = 0.0 # flag = ALIVE=1 or DEAD=0 \nReflBin = [None] * Nbinsp1 #bin to store weights of escaped photos for reflectivity\nepi_albedo = epi_mus\/(epi_mus + epi_mua) # albedo of tissue\nderm_albedo = derm_mus\/(derm_mus + derm_mua) # albedo of tissue\nNphotons = 1000 # number of photons in simulation \nNR = Nbins # number of radial positions \nradial_size = 2.5 # maximum radial size \nr = 0.0 # radial position \ndr = radial_size\/NR; # cm, radial bin size \nir = 0 # index to radial position \nshellvolume = 0.0 # volume of shell at radial position r \nCNT = 0.0 # total count of photon weight summed over all bins \nrnd = 0.0 # assigned random value 0-1 \nu = 0.0\ntemp = 0.0 # dummy variables<\/code><\/pre>\n\n\n\n<\/p>\n\n\n\n
And the last thing before we begin looping; some initialization:<\/p>\n\n\n\n
# Inits\nrandom.seed(0)\nRandomNum = random.random()\nfor i in range(NR+1):\n ReflBin[i] = 0<\/code><\/pre>\n\n\n\nNow we can begin our iteration code. We\u2019ll start with a simple while loop for number of photons:<\/p>\n\n\n\n
while True:\n i_photon = i_photon + 1\n if i_photon >= Nphotons:\n break<\/code><\/pre>\n\n\n\nWithin <\/em>this loop, we need to do some further initialization for this loop:<\/p>\n\n\n\nW = 1.0
photon_status = ALIVE
x= 0
y = 0
z = 0
#Randomly set photon trajectory to yield an isotropic source.
costheta = 2.0 * random.random() – 1.0
sintheta = math.sqrt(1.0 – costheta*costheta)
psi = 2.0 * PI * random.random()
ux = sintheta * math.cos(psi)
uy = sintheta * math.sin(psi)
uz = (abs(costheta)) # on the first step we want to head down, into the tissue, so > 0
# Propagate one photon until it dies as determined by ROULETTE.
# or if it reaches the surface again
it = 0
max_iterations = 100000 # to help avoid infinite loops in case we do something wrong
# we’ll hit epidermis first, so set mua\/mus to those scattering\/absorption values
mua = epi_mua
mus = epi_mus
albedo = epi_albedo
<\/p>\n\n\n\n
W = 1.0
\nphoton_status = ALIVE<\/code>\n<\/code><\/p>\n\n\n\nx= 0
\ny = 0
\nz = 0<\/code><\/p>\n\n\n\n\n#Randomly set photon trajectory to yield an isotropic source.
\ncostheta = 2.0 * random.random() - 1.0<\/p>\n
sintheta = math.sqrt(1.0 - costheta*costheta)
\npsi = 2.0 * PI * random.random()
\nux = sintheta * math.cos(psi)
\nuy = sintheta * math.sin(psi)
\nuz = (abs(costheta)) # on the first step we want to head down, into the tissue, so > 0<\/p>\n
# Propagate one photon until it dies as determined by ROULETTE.
\n# or if it reaches the surface again
\nit = 0
\nmax_iterations = 100000 # to help avoid infinite loops in case we do something wrong<\/p>\n
# we'll hit epidermis first, so set mua\/mus to those scattering\/absorption values
\nmua = epi_mua
\nmus = epi_mus
\nalbedo = epi_albedo<\/p>\n<\/code>\n\n\n\n
<\/p>\n\n\n\n
W = 1.0
\nphoton_status = ALIVE<\/code><\/p>\n\n\n\nW = 1.0
\nphoton_status = ALIVE<\/code>
\n<\/p>\n\n\n\nx= 0
\ny = 0
\nz = 0<\/p>\n\n\n\n
#Randomly set photon trajectory to yield an isotropic source.
\ncostheta = 2.0 * random.random() – 1.0<\/p>\n\n\n\n
sintheta = math.sqrt(1.0 – costheta*costheta)
\npsi = 2.0 * PI * random.random()
\nux = sintheta * math.cos(psi)
\nuy = sintheta * math.sin(psi)
\nuz = (abs(costheta)) # on the first step we want to head down, into the tissue, so > 0<\/p>\n\n\n\n
# Propagate one photon until it dies as determined by ROULETTE.
\n# or if it reaches the surface again
\nit = 0
\nmax_iterations = 100000 # to help avoid infinite loops in case we do something wrong<\/p>\n\n\n\n
# we’ll hit epidermis first, so set mua\/mus to those scattering\/absorption values
\nmua = epi_mua
\nmus = epi_mus
\nalbedo = epi_albedo<\/p>\n\n\n\n
This will so far spawn a photon and do nothing. We now need to implement that flow chart from the Theory part; launch, hop, drop, spin, terminate.<\/p>\n\n\n\n
Let\u2019s make a new loop within the existing one. The first thing we have to do is add a break condition \u2013 for now, we\u2019ll put in a max iteration, but we\u2019ll also add in our terminate functionality.<\/p>\n\n\n\n
while True: \n it = it + 1\n # Check Roulette\n if (W < THRESHOLD):\n if (random.random() <= CHANCE):\n W = W \/ CHANCE\n else:\n photon_status = DEAD\n if photon_status is DEAD:\n break\n if it > max_iterations:\n break<\/code><\/pre>\n\n\n\nwhile True:
\n it = it + 1
\n # Check Roulette
\n if (W < THRESHOLD):
\n if (random.random() <= CHANCE):
\n W = W \/ CHANCE
\n else:
\n photon_status = DEAD
\n if photon_status is DEAD:
\n break
\n if it > max_iterations:
\n break<\/code><\/p>\n\n\n\nwhile True:
\n it = it + 1
\n # Check Roulette
\n if (W < THRESHOLD):
\n if (random.random() <= CHANCE):
\n W = W \/ CHANCE
\n else:
\n photon_status = DEAD
\n if photon_status is DEAD:
\n break
\n if it > max_iterations:
\n break<\/code><\/p>\n\n\n\nwhile True:
\n it = it + 1
\n # Check Roulette
\n if (W < THRESHOLD):
\n if (random.random() <= CHANCE):
\n W = W \/ CHANCE
\n else:
\n photon_status = DEAD
\n if photon_status is DEAD:
\n break
\n if it > max_iterations:
\n break<\/code><\/p>\n\n\n\nAbove the \u2018# Check Roulette<\/em>\u2018 line, we\u2019ll flesh out our hop code:<\/p>\n\n\n\n# Check Roulette<\/em><\/p>\n\n\n\n# Check Roulette<\/em><\/p>\n\n\n\n# Check Roulette<\/em><\/p>\n\n\n\nrnd = random.random()\nwhile rnd <= 0.0: # make sure it is > 0.0\n rnd = random.random()\ns = -math.log(rnd)\/(mua + mus)\nx = x + (s * ux)\ny = y + (s * uy)\nz = z + (s * uz)<\/code><\/pre>\n\n\n\nrnd = random.random()
\nwhile rnd <= 0.0: # make sure it is > 0.0
\n rnd = random.random()
\ns = -math.log(rnd)\/(mua + mus)
\nx = x + (s * ux)
\ny = y + (s * uy)
\nz = z + (s * uz)<\/code><\/p>\n\n\n\nrnd = random.random()
\nwhile rnd <= 0.0: # make sure it is > 0.0
\n rnd = random.random()
\ns = -math.log(rnd)\/(mua + mus)
\nx = x + (s * ux)
\ny = y + (s * uy)
\nz = z + (s * uz)<\/code><\/p>\n\n\n\nrnd = random.random()
\nwhile rnd <= 0.0: # make sure it is > 0.0
\n rnd = random.random()
\ns = -math.log(rnd)\/(mua + mus)
\nx = x + (s * ux)
\ny = y + (s * uy)
\nz = z + (s * uz)<\/code><\/p>\n\n\n\nOn our first \u2018hop\u2019, the code is the \u2018launch\u2019 step. Immediately after this, we do a boundary check to see if the photon is escaping the skin, moving back to the boundary as needed:<\/p>\n\n\n\n
if uz < 0:\n # calculate partial step to reach boundary surface\n s1 = abs(z\/uz)\n # move back\n x = x - (s * ux)\n y = y - (s * uy)\n z = z - (s * uz)\n # take partial step\n x = x + (s1 * ux)\n y = y + (s1 * uy)\n z = z + (s1 * uz)<\/code><\/pre>\n\n\n\nif uz < 0:
\n # calculate partial step to reach boundary surface
\n s1 = abs(z\/uz)
\n # move back
\n x = x - (s * ux)
\n y = y - (s * uy)
\n z = z - (s * uz)
\n # take partial step
\n x = x + (s1 * ux)
\n y = y + (s1 * uy)
\n z = z + (s1 * uz)<\/code><\/p>\n\n\n\nif uz < 0:
\n # calculate partial step to reach boundary surface
\n s1 = abs(z\/uz)
\n # move back
\n x = x - (s * ux)
\n y = y - (s * uy)
\n z = z - (s * uz)
\n # take partial step
\n x = x + (s1 * ux)
\n y = y + (s1 * uy)
\n z = z + (s1 * uz)<\/code><\/p>\n\n\n\nif uz < 0:
\n # calculate partial step to reach boundary surface
\n s1 = abs(z\/uz)
\n # move back
\n x = x - (s * ux)
\n y = y - (s * uy)
\n z = z - (s * uz)
\n # take partial step
\n x = x + (s1 * ux)
\n y = y + (s1 * uy)
\n z = z + (s1 * uz)<\/code><\/p>\n\n\n\nTo finish of this if-statement, we need to find out how much of our photon energy escaped and how much was reflected internally. To do this, we call our Fresnel function, and bounce the photon back:<\/p>\n\n\n\n
# photon is now at the surface boundary, figure out how much escaped and how much was reflected\ninternal_reflectance = self.RFresnel(1.0,nt, -uz )\n\n#Add weighted reflectance of escpaed photon to reflectance bin\nexternal_reflectance = 1 - internal_reflectance\nr = math.sqrt(x*x + y*y)\nir = (r\/dr)\nif ir >= NR:\n ir = NR \nReflBin[int(ir)] = ReflBin[int(ir)] + (W * external_reflectance)\n\t\t\t\t\t\t\n# Bounce the photon back into the skin\nW = internal_reflectance * W\nuz = -uz\nx = (s-s1) * ux\ny = (s-s1) * uy\nz = (s-s1) * uz<\/code><\/pre>\n\n\n\n# photon is now at the surface boundary, figure out how much escaped and how much was reflected
\ninternal_reflectance = self.RFresnel(1.0,nt, -uz )<\/code>\n<\/code><\/p>\n\n\n\n#Add weighted reflectance of escpaed photon to reflectance bin
\nexternal_reflectance = 1 - internal_reflectance
\nr = math.sqrt(x*x + y*y)
\nir = (r\/dr)
\nif ir >= NR:
\n ir = NR
\nReflBin[int(ir)] = ReflBin[int(ir)] + (W * external_reflectance)<\/code><\/p>\n\n\n\n\n<\/code>\n\n\n\n<\/p>\n\n\n\n
# Bounce the photon back into the skin
\nW = internal_reflectance * W
\nuz = -uz
\nx = (s-s1) * ux
\ny = (s-s1) * uy
\nz = (s-s1) * uz<\/code><\/p>\n\n\n\n# photon is now at the surface boundary, figure out how much escaped and how much was reflected
\ninternal_reflectance = self.RFresnel(1.0,nt, -uz )<\/code><\/p>\n\n\n\n# photon is now at the surface boundary, figure out how much escaped and how much was reflected
\ninternal_reflectance = self.RFresnel(1.0,nt, -uz )<\/code>
\n<\/p>\n\n\n\n#Add weighted reflectance of escpaed photon to reflectance bin
\nexternal_reflectance = 1 – internal_reflectance
\nr = math.sqrt(x*x + y*y)
\nir = (r\/dr)
\nif ir >= NR:
\n ir = NR
\nReflBin[int(ir)] = ReflBin[int(ir)] + (W * external_reflectance)<\/p>\n\n\n\n
# Bounce the photon back into the skin
\nW = internal_reflectance * W
\nuz = -uz
\nx = (s-s1) * ux
\ny = (s-s1) * uy
\nz = (s-s1) * uz<\/p>\n\n\n\n
Close of the if-statement. Rather than do a full boundary check between the skin layers, we\u2019re simply going to adjust our albedo and scattering parameters. It makes the code a bit easier to maintain, but if you want more of a challenge, you can also put reflectance between the skin layers.<\/p>\n\n\n\n
# check if we have passed into the second layer, or the first\nif z <= epidermis_thickness:\n mua = epi_mua\n mus = epi_mus\n albedo = epi_albedo\nelse:\n mua = derm_mua\n mus = derm_mus\n albedo = derm_albedo<\/code><\/pre>\n\n\n\n# check if we have passed into the second layer, or the first
\nif z <= epidermis_thickness:
\n mua = epi_mua
\n mus = epi_mus
\n albedo = epi_albedo
\nelse:
\n mua = derm_mua
\n mus = derm_mus
\n albedo = derm_albedo<\/code><\/p>\n\n\n\n# check if we have passed into the second layer, or the first
\nif z <= epidermis_thickness:
\n mua = epi_mua
\n mus = epi_mus
\n albedo = epi_albedo
\nelse:
\n mua = derm_mua
\n mus = derm_mus
\n albedo = derm_albedo<\/code><\/p>\n\n\n\n# check if we have passed into the second layer, or the first
\nif z <= epidermis_thickness:
\n mua = epi_mua
\n mus = epi_mus
\n albedo = epi_albedo
\nelse:
\n mua = derm_mua
\n mus = derm_mus
\n albedo = derm_albedo<\/code><\/p>\n\n\n\nDrop stage now \u2013 fairly simple. We just reduce the photon weight by the albedo.<\/p>\n\n\n\n
absorb = W*(1 - albedo)\nW = W - absorb<\/code><\/pre>\n\n\n\nabsorb = W*(1 - albedo)
\nW = W - absorb<\/code><\/p>\n\n\n\nabsorb = W*(1 - albedo)
\nW = W - absorb<\/code><\/p>\n\n\n\nabsorb = W*(1 - albedo)
\nW = W - absorb<\/code>
\n<\/p>\n\n\n\nStraight into the Spin\u2026 :<\/p>\n\n\n\n
''' SPIN \n Scatter photon into new trajectory defined by theta and psi.\n Theta is specified by cos(theta), which is determined \n based on the Henyey-Greenstein scattering function.\n Convert theta and psi into cosines ux, uy, uz. \n'''\n# Sample for costheta \nrnd = random.random()\nif (g == 0.0):\n costheta = 2.0*rnd - 1.0\nelse:\n temp = (1.0 - g*g)\/(1.0 - g + 2*g*rnd)\n costheta = (1.0 + g*g - temp*temp)\/(2.0*g)\nsintheta = math.sqrt(1.0 - costheta*costheta) \n\n# Sample psi. \npsi = 2.0*PI*random.random()\ncospsi = math.cos(psi)\nif (psi < PI):\n sinpsi = math.sqrt(1.0 - cospsi*cospsi)\nelse:\n sinpsi = -math.sqrt(1.0 - cospsi*cospsi)\n\n# New trajectory. \nif (1 - abs(uz) <= ONE_MINUS_COSZERO) : # close to perpendicular. \n uxx = sintheta * cospsi\n uyy = sintheta * sinpsi\n uzz = costheta * self.SIGN(uz) # SIGN() is faster than division. \nelse: # usually use this option \n temp = math.sqrt(1.0 - uz * uz)\n uxx = sintheta * (ux * uz * cospsi - uy * sinpsi) \/ temp + ux * costheta\n uyy = sintheta * (uy * uz * cospsi + ux * sinpsi) \/ temp + uy * costheta\n uzz = -sintheta * cospsi * temp + uz * costheta\n\n# Update trajectory \nux = uxx\nuy = uyy\nuz = uzz<\/code><\/pre>\n\n\n\n''' SPIN
\n Scatter photon into new trajectory defined by theta and psi.
\n Theta is specified by cos(theta), which is determined
\n based on the Henyey-Greenstein scattering function.
\n Convert theta and psi into cosines ux, uy, uz.
\n'''
\n# Sample for costheta
\nrnd = random.random()
\nif (g == 0.0):
\n costheta = 2.0*rnd - 1.0
\nelse:
\n temp = (1.0 - g*g)\/(1.0 - g + 2*g*rnd)
\n costheta = (1.0 + g*g - temp*temp)\/(2.0*g)
\nsintheta = math.sqrt(1.0 - costheta*costheta) <\/code>\n<\/code><\/p>\n\n\n\n# Sample psi.
\npsi = 2.0*PI*random.random()
\ncospsi = math.cos(psi)
\nif (psi < PI):
\n sinpsi = math.sqrt(1.0 - cospsi*cospsi)
\nelse:
\n sinpsi = -math.sqrt(1.0 - cospsi*cospsi)<\/code><\/p>\n\n\n\n\n# New trajectory.
\nif (1 - abs(uz) <= ONE_MINUS_COSZERO) : # close to perpendicular.
\n uxx = sintheta * cospsi
\n uyy = sintheta * sinpsi
\n uzz = costheta * self.SIGN(uz) # SIGN() is faster than division.
\nelse: # usually use this option
\n temp = math.sqrt(1.0 - uz * uz)
\n uxx = sintheta * (ux * uz * cospsi - uy * sinpsi) \/ temp + ux * costheta
\n uyy = sintheta * (uy * uz * cospsi + ux * sinpsi) \/ temp + uy * costheta
\n uzz = -sintheta * cospsi * temp + uz * costheta<\/p>\n<\/code>\n\n\n\n
<\/p>\n\n\n\n
# Update trajectory
\nux = uxx
\nuy = uyy
\nuz = uzz<\/code><\/p>\n\n\n\n''' SPIN
\n Scatter photon into new trajectory defined by theta and psi.
\n Theta is specified by cos(theta), which is determined
\n based on the Henyey-Greenstein scattering function.
\n Convert theta and psi into cosines ux, uy, uz.
\n'''
\n# Sample for costheta
\nrnd = random.random()
\nif (g == 0.0):
\n costheta = 2.0*rnd - 1.0
\nelse:
\n temp = (1.0 - g*g)\/(1.0 - g + 2*g*rnd)
\n costheta = (1.0 + g*g - temp*temp)\/(2.0*g)
\nsintheta = math.sqrt(1.0 - costheta*costheta) <\/code><\/p>\n\n\n\n''' SPIN
\n Scatter photon into new trajectory defined by theta and psi.
\n Theta is specified by cos(theta), which is determined
\n based on the Henyey-Greenstein scattering function.
\n Convert theta and psi into cosines ux, uy, uz.
\n'''
\n# Sample for costheta
\nrnd = random.random()
\nif (g == 0.0):
\n costheta = 2.0*rnd - 1.0
\nelse:
\n temp = (1.0 - g*g)\/(1.0 - g + 2*g*rnd)
\n costheta = (1.0 + g*g - temp*temp)\/(2.0*g)
\nsintheta = math.sqrt(1.0 - costheta*costheta) <\/code><\/p>\n\n\n\n# Sample psi.
\npsi = 2.0*PI*random.random()
\ncospsi = math.cos(psi)
\nif (psi < PI):
\n sinpsi = math.sqrt(1.0 – cospsi*cospsi)
\nelse:
\n sinpsi = -math.sqrt(1.0 – cospsi*cospsi)<\/p>\n\n\n\n
# New trajectory.
\nif (1 – abs(uz) <= ONE_MINUS_COSZERO) : # close to perpendicular.
\n uxx = sintheta * cospsi
\n uyy = sintheta * sinpsi
\n uzz = costheta * self.SIGN(uz) # SIGN() is faster than division.
\nelse: # usually use this option
\n temp = math.sqrt(1.0 – uz * uz)
\n uxx = sintheta * (ux * uz * cospsi – uy * sinpsi) \/ temp + ux * costheta
\n uyy = sintheta * (uy * uz * cospsi + ux * sinpsi) \/ temp + uy * costheta
\n uzz = -sintheta * cospsi * temp + uz * costheta<\/p>\n\n\n\n
# Update trajectory
\nux = uxx
\nuy = uyy
\nuz = uzz<\/p>\n\n\n\n
The terminate code we wrote towards the beginning should sit below all of the above, being the final stage.<\/p>\n\n\n\n
Calculating Reflectance<\/h2>\n\n\n\n
The code should run, but nothing will happen yet. We need to sum up the reflections, as described in the previous chapter. We will add this beneath our nested while loops at the end of our Monte Carlo function:<\/p>\n\n\n\n
total_reflection = 0.0\n for each in range(NR+1):\n total_reflection = total_reflection + ReflBin[each]\/Nphotons\n return total_reflection<\/code><\/pre>\n\n\n\ntotal_reflection = 0.0
\n for each in range(NR+1):
\n total_reflection = total_reflection + ReflBin[each]\/Nphotons
\n return total_reflection<\/code><\/p>\n\n\n\ntotal_reflection = 0.0
\n for each in range(NR+1):
\n total_reflection = total_reflection + ReflBin[each]\/Nphotons
\n return total_reflection<\/code><\/p>\n\n\n\ntotal_reflection = 0.0
\n for each in range(NR+1):
\n total_reflection = total_reflection + ReflBin[each]\/Nphotons
\n return total_reflection<\/code><\/p>\n\n\n\nNow we need up update our GetReflectanceValues<\/em>() function so that it can manage the data returned. For debugging purposes, we\u2019ll put the reflectance values on the end of our CSV debug code, like so:<\/p>\n\n\n\nGetReflectanceValues<\/em>
\n
\nGetReflectanceValues<\/em><\/p>\n\n\n\n[....]\n# Scattering coefficients\nscattering_epidermis = pow(14.74 * nm, -0.22) + 2.2 * pow(10,11) * pow(nm, -4)\nscattering_dermis = scattering_epidermis * 0.5\n \nreflectance = self.MonteCarlo(epidermis, scattering_epidermis, dermis, scattering_dermis, nm)\n \nCSV_Out = CSV_Out + f\"{nm},{SAV_eumelanin_L},{SAV_pheomelanin_L},{baselineSkinAbsorption_L},{epidermis},{scattering_epidermis},{scattering_dermis},{dermis},{reflectance}\\n\"\n return CSV_Out\n[...]<\/code><\/pre>\n\n\n\n[....]
\n# Scattering coefficients
\nscattering_epidermis = pow(14.74 * nm, -0.22) + 2.2 * pow(10,11) * pow(nm, -4)
\nscattering_dermis = scattering_epidermis * 0.5<\/code>\n<\/code><\/p>\n\n\n\nreflectance = self.MonteCarlo(epidermis, scattering_epidermis, dermis, scattering_dermis, nm)<\/code><\/p>\n\n\n\n\n<\/code>\n\n\n\n<\/p>\n\n\n\n
CSV_Out = CSV_Out + f\"{nm},{SAV_eumelanin_L},{SAV_pheomelanin_L},{baselineSkinAbsorption_L},{epidermis},{scattering_epidermis},{scattering_dermis},{dermis},{reflectance}\\n\"
\n return CSV_Out
\n[...]<\/code><\/p>\n\n\n\n[....]
\n# Scattering coefficients
\nscattering_epidermis = pow(14.74 * nm, -0.22) + 2.2 * pow(10,11) * pow(nm, -4)
\nscattering_dermis = scattering_epidermis * 0.5<\/code><\/p>\n\n\n\n[....]
\n# Scattering coefficients
\nscattering_epidermis = pow(14.74 * nm, -0.22) + 2.2 * pow(10,11) * pow(nm, -4)
\nscattering_dermis = scattering_epidermis * 0.5<\/code><\/p>\n\n\n\nreflectance = self.MonteCarlo(epidermis, scattering_epidermis, dermis, scattering_dermis, nm)<\/p>\n\n\n\n
CSV_Out = CSV_Out + f”{nm},{SAV_eumelanin_L},{SAV_pheomelanin_L},{baselineSkinAbsorption_L},{epidermis},{scattering_epidermis},{scattering_dermis},{dermis},{reflectance}\\n”
\n return CSV_Out
\n[…]<\/p>\n\n\n\n
Also don\u2019t forget to update the headers in Generate()<\/em>\u2026<\/p>\n\n\n\nGenerate()<\/em><\/p>\n\n\n\nGenerate()<\/em><\/p>\n\n\n\nGenerate()<\/em><\/p>\n\n\n\ndef Generate(self):\n CSV_Out = \"Wavelength,Eumelanin,Pheomelanin,Baseline,Epidermis,ScatteringEpi,ScatteringDerm,Dermis,Reflectance\\n\"\n<\/code><\/pre>\n\n\n\ndef Generate(self):
\n CSV_Out = \"Wavelength,Eumelanin,Pheomelanin,Baseline,Epidermis,ScatteringEpi,ScatteringDerm,Dermis,Reflectance\\n\"
\n<\/code><\/p>\n\n\n\ndef Generate(self):
\n CSV_Out = \"Wavelength,Eumelanin,Pheomelanin,Baseline,Epidermis,ScatteringEpi,ScatteringDerm,Dermis,Reflectance\\n\"
\n<\/code><\/p>\n\n\n\ndef Generate(self):
\n CSV_Out = \"Wavelength,Eumelanin,Pheomelanin,Baseline,Epidermis,ScatteringEpi,ScatteringDerm,Dermis,Reflectance\\n\"
\n<\/code><\/p>\n\n\n\nDebugging<\/h2>\n\n\n\n
Give the code a run and import your CSV into Excel again.This may take some time<\/em><\/span><\/strong>, it takes a good 15 mins on a Core i9-9900k. I haven\u2019t implemented multi-threading for this tutorial as it complicates what is already a heavy topic. If you\u2019re confident modifying this to be multi-threaded, feel free to do so. <\/p>\n\n\n\nThis may take some time<\/em><\/span><\/strong>
\nThis may take some time<\/em><\/span>
\nThis may take some time<\/em><\/p>\n\n\n\nThis may take some time<\/em><\/span><\/strong>
\nThis may take some time<\/em><\/span>
\nThis may take some time<\/em><\/p>\n\n\n\nThis may take some time<\/em><\/span><\/strong>
\nThis may take some time<\/em><\/span>
\nThis may take some time<\/em>
\n
\nThis may take some time<\/em><\/span>
\nThis may take some time<\/em>
\n
\nThis may take some time<\/em><\/p>\n\n\n\nThe final code:<\/p>\n\n\n\n
import random\nimport math\n\nclass SkinAbsorption:\n\n O2Hb = {}\n Hb = {}\n \n def __init__(self):\n print(\"init SkinAbsorption\")\n HbFile = open('hb.csv', \"r\")\n O2HbFile = open('O2Hb.csv', \"r\")\n\n HbLines = HbFile.readlines()\n for line in HbLines:\n splitLine = line.split(\",\")\n self.Hb[int(splitLine[0])] = float(splitLine[1].rstrip(\"\\n\"))\n\n O2HbLines = O2HbFile.readlines()\n for line in O2HbLines:\n splitLine = line.split(\",\")\n self.O2Hb[int(splitLine[0])] = float(splitLine[1].rstrip(\"\\n\"))\n \n HbFile.close()\n O2HbFile.close()\n\n def Generate(self):\n CSV_Out = \"Wavelength,Cm, Ch, Bm,Eumelanin,Pheomelanin,Baseline,Epidermis,ScatteringEpi,ScatteringDerm,Dermis,Reflectance\\n\"\n groupByBlend = []\n for Bm in [0.01,0.5,0.99]:\n groupByMelanin = []\n for Cm in [0.002,0.0135,0.0425,0.1,0.185,0.32,0.5]:\n groupByHemoglobin = []\n for Ch in [0.003,0.02,0.07,0.16,0.32]:\n values = self.GetReflectanceValues( Cm, Ch, Bm )\n CSV_Out = CSV_Out + values\n groupByHemoglobin.append(values)\n groupByMelanin.append(groupByHemoglobin)\n groupByBlend.append(groupByMelanin)\n CSV_Out_File = open(\"AbsorptionValues.csv\",\"w+\")\n CSV_Out_File.write(CSV_Out)\n CSV_Out_File.close() \n return groupByBlend\n \n def GetReflectanceValues(self, Cm, Ch, Bm):\n CSV_Out = \"\"\n wavelengths = range(380,790,10)\n for nm in wavelengths:\n \n # First layer absorption - Epidermis\n SAV_eumelanin_L = 6.6 * pow(10,10) * pow(nm,-3.33)\n SAV_pheomelanin_L = 2.9 * pow(10,14) * pow(nm,-4.75)\n epidermal_hemoglobin_fraction = Ch * 0.25\n \n # baseline - used in both layers\n baselineSkinAbsorption_L = 0.0244 + 8.54 * pow(10,-(nm-220)\/123)\n\n # Epidermis Absorption Coefficient: \n epidermis = Cm * ((Bm * SAV_eumelanin_L) +((1 - Bm ) * SAV_pheomelanin_L)) + ((1 - epidermal_hemoglobin_fraction) * baselineSkinAbsorption_L)\n \n # Second layer absorption - Dermis\n gammaOxy = self.O2Hb[int(nm)]\n gammaDeoxy = self.Hb[int(nm)]\n A = 0.75 * gammaOxy + (1 - 0.75) * gammaDeoxy \n dermis = Ch * (A + (( 1 - Ch) * baselineSkinAbsorption_L))\n \n # Scattering coefficients\n scattering_epidermis = pow(14.74 * nm, -0.22) + 2.2 * pow(10,11) * pow(nm, -4)\n scattering_dermis = scattering_epidermis * 0.5\n \n reflectance = self.MonteCarlo(epidermis, scattering_epidermis, dermis, scattering_dermis, nm)\n print(f\"Done on {nm} >> ({Cm}, {Ch}, {Bm}) = {reflectance}\")\n \n CSV_Out = CSV_Out + f\"{nm},{Cm},{Ch},{Bm},{SAV_eumelanin_L},{SAV_pheomelanin_L},{baselineSkinAbsorption_L},{epidermis},{scattering_epidermis},{scattering_dermis},{dermis},{reflectance}\\n\"\n return CSV_Out\n\n def MonteCarlo (self, epi_mua, epi_mus, derm_mua, derm_mus, nm):\n # These are our Monte Carlo Light Transport Variables that don't change\n Nbins = 1000\n Nbinsp1 = 1001\n PI = 3.1415926\n LIGHTSPEED = 2.997925 * pow(10,10)\n ALIVE = 1\n DEAD = 0\n THRESHOLD = 0.01\n CHANCE = 0.1\n COS90D = 1 * pow(10,-6)\n ONE_MINUS_COSZERO = 1 * pow(10,-12)\n COSZERO = 1.0 - 1.0e-12 # cosine of about 1e-6 rad\n g = 0.9\n nt = 1.33 # Index of refraction\n epidermis_thickness = 0.25\n\n x = 0.0\n y = 0.0\n z = 0.0 # photon position\n\n ux = 0.0\n uy = 0.0\n uz = 0.0 # photon trajectory as cosines\n\n uxx = 0.0\n uyy = 0.0\n uzz = 0.0 # temporary values used during SPIN\n\n s = 0.0 # step sizes. s = -log(RND)\/mus [cm] \n costheta = 0.0 # cos(theta) \n sintheta = 0.0 # sin(theta) \n cospsi = 0.0 # cos(psi) \n sinpsi = 0.0 # sin(psi) \n psi = 0.0 # azimuthal angle \n i_photon = 0.0 # current photon\n W = 0.0 # photon weight \n absorb = 0.0 # weighted deposited in a step due to absorption \n photon_status = 0.0 # flag = ALIVE=1 or DEAD=0 \n ReflBin = [None] * Nbinsp1 #bin to store weights of escaped photos for reflectivity\n epi_albedo = epi_mus\/(epi_mus + epi_mua) # albedo of tissue\n derm_albedo = derm_mus\/(derm_mus + derm_mua) # albedo of tissue\n Nphotons = 1000 # number of photons in simulation \n NR = Nbins # number of radial positions \n radial_size = 2.5 # maximum radial size \n r = 0.0 # radial position \n dr = radial_size\/NR; # cm, radial bin size \n ir = 0 # index to radial position \n shellvolume = 0.0 # volume of shell at radial position r \n CNT = 0.0 # total count of photon weight summed over all bins \n rnd = 0.0 # assigned random value 0-1 \n u = 0.0\n temp = 0.0 # dummy variables\n\n # Inits\n random.seed(0)\n RandomNum = random.random()\n for i in range(NR+1):\n ReflBin[i] = 0\n\n while True:\n i_photon = i_photon + 1\n\n W = 1.0\n photon_status = ALIVE\n\n x= 0\n y = 0\n z = 0\n\n #Randomly set photon trajectory to yield an isotropic source.\n costheta = 2.0 * random.random() - 1.0\n\n sintheta = math.sqrt(1.0 - costheta*costheta)\n psi = 2.0 * PI * random.random()\n ux = sintheta * math.cos(psi)\n uy = sintheta * math.sin(psi)\n uz = (abs(costheta)) # on the first step we want to head down, into the tissue, so > 0\n\n # Propagate one photon until it dies as determined by ROULETTE.\n # or if it reaches the surface again\n it = 0\n max_iterations = 100000 # to help avoid infinite loops in case we do something wrong\n\n # we'll hit epidermis first, so set mua\/mus to those scattering\/absorption values\n mua = epi_mua\n mus = epi_mus\n albedo = epi_albedo\n while True: \n it = it + 1\n rnd = random.random()\n while rnd <= 0.0: # make sure it is > 0.0\n rnd = random.random()\n s = -math.log(rnd)\/(mua + mus)\n x = x + (s * ux)\n y = y + (s * uy)\n z = z + (s * uz)\n\n if uz < 0:\n # calculate partial step to reach boundary surface\n s1 = abs(z\/uz)\n # move back\n x = x - (s * ux)\n y = y - (s * uy)\n z = z - (s * uz)\n # take partial step\n x = x + (s1 * ux)\n y = y + (s1 * uy)\n z = z + (s1 * uz)\n\n # photon is now at the surface boundary, figure out how much escaped and how much was reflected\n internal_reflectance = self.RFresnel(1.0,nt, -uz )\n\n #Add weighted reflectance of escpaed photon to reflectance bin\n external_reflectance = 1 - internal_reflectance\n r = math.sqrt(x*x + y*y)\n ir = (r\/dr)\n if ir >= NR:\n ir = NR \n ReflBin[int(ir)] = ReflBin[int(ir)] + (W * external_reflectance)\n \n # Bounce the photon back into the skin\n W = internal_reflectance * W\n uz = -uz\n x = (s-s1) * ux\n y = (s-s1) * uy\n z = (s-s1) * uz\n\n # check if we have passed into the second layer, or the first\n if z <= epidermis_thickness:\n mua = epi_mua\n mus = epi_mus\n albedo = epi_albedo\n else:\n mua = derm_mua\n mus = derm_mus\n albedo = derm_albedo\n\n ''' DROP '''\n absorb = W*(1 - albedo)\n W = W - absorb\n\n ''' SPIN '''\n # Sample for costheta \n rnd = random.random()\n if (g == 0.0):\n costheta = 2.0*rnd - 1.0\n else:\n temp = (1.0 - g*g)\/(1.0 - g + 2*g*rnd)\n costheta = (1.0 + g*g - temp*temp)\/(2.0*g)\n sintheta = math.sqrt(1.0 - costheta*costheta) \n\n # Sample psi. \n psi = 2.0*PI*random.random()\n cospsi = math.cos(psi)\n if (psi < PI):\n sinpsi = math.sqrt(1.0 - cospsi*cospsi)\n else:\n sinpsi = -math.sqrt(1.0 - cospsi*cospsi)\n\n # New trajectory. \n if (1 - abs(uz) <= ONE_MINUS_COSZERO) : # close to perpendicular. \n uxx = sintheta * cospsi\n uyy = sintheta * sinpsi\n uzz = costheta * self.SIGN(uz) # SIGN() is faster than division. \n else: # usually use this option \n temp = math.sqrt(1.0 - uz * uz)\n uxx = sintheta * (ux * uz * cospsi - uy * sinpsi) \/ temp + ux * costheta\n uyy = sintheta * (uy * uz * cospsi + ux * sinpsi) \/ temp + uy * costheta\n uzz = -sintheta * cospsi * temp + uz * costheta\n\n # Update trajectory \n ux = uxx\n uy = uyy\n uz = uzz\n\n \n # Check Roulette\n if (W < THRESHOLD):\n if (random.random() <= CHANCE):\n W = W \/ CHANCE\n else:\n photon_status = DEAD\n if photon_status is DEAD:\n break\n if it > max_iterations:\n break\n \n if i_photon >= Nphotons:\n break\n total_reflection = 0.0\n for each in range(NR+1):\n total_reflection = total_reflection + ReflBin[each]\/Nphotons\n return total_reflection\n\n def RFresnel(self, n1, n2, cosT1):\n r = 0.0\n cosT2 = 0.0\n COSZERO = 1.0 - 1.0e-12\n COS90D = 1 * pow(10,-6)\n if n1 == n2: #matched boundary\n r = 0.0\n cosT2 = cosT1\n elif cosT1 > COSZERO: # normal incident\n cosT2 = ca1\n r = (n2-n1)\/(n2+n1)\n r *= r\n elif cosT1 < COS90D: # very slant\n cosT2 = 0.0\n r = 1.0\n else: #general\n sinT1 = math.sqrt(1 - cosT1*cosT1)\n sinT2 = n1 * sinT1\/n2\n \n if sinT2 >= 1.0:\n r = 1.0\n cosT2 = 0.0\n else:\n cosT2 = math.sqrt(1 - sinT2 * sinT2)\n cosAP = cosT1*cosT2 - sinT1*sinT2\n cosAM = cosT1*cosT2 + sinT1*sinT2\n sinAP = sinT1*cosT2 + cosT1*sinT2\n sinAM = sinT1*cosT2 - cosT1*sinT2\n r = 0.5 * sinAM * sinAM*(cosAM*cosAM+cosAP*cosAP)\/(sinAP*sinAP*cosAM*cosAM)\n return r\n\n def SIGN(self, x):\n if x >=0:\n return 1\n else:\n return 0 \n\nskinAbsorption = SkinAbsorption()\nreflectanceValues = skinAbsorption.Generate()<\/code><\/pre>\n\n\n\nimport random
\nimport math<\/code>\n<\/code><\/p>\n\n\n\nclass SkinAbsorption:<\/code><\/p>\n\n\n\n\n O2Hb = {}
\n Hb = {}<\/p>\n
def __init__(self):
\n print(\"init SkinAbsorption\")
\n HbFile = open('hb.csv', \"r\")
\n O2HbFile = open('O2Hb.csv', \"r\")<\/p>\n
HbLines = HbFile.readlines()
\n for line in HbLines:
\n splitLine = line.split(\",\")
\n self.Hb[int(splitLine[0])] = float(splitLine[1].rstrip(\"\\n\"))<\/p>\n
O2HbLines = O2HbFile.readlines()
\n for line in O2HbLines:
\n splitLine = line.split(\",\")
\n self.O2Hb[int(splitLine[0])] = float(splitLine[1].rstrip(\"\\n\"))<\/p>\n
HbFile.close()
\n O2HbFile.close()<\/p>\n
def Generate(self):
\n CSV_Out = \"Wavelength,Cm, Ch, Bm,Eumelanin,Pheomelanin,Baseline,Epidermis,ScatteringEpi,ScatteringDerm,Dermis,Reflectance\\n\"
\n groupByBlend = []
\n for Bm in [0.01,0.5,0.99]:
\n groupByMelanin = []
\n for Cm in [0.002,0.0135,0.0425,0.1,0.185,0.32,0.5]:
\n groupByHemoglobin = []
\n for Ch in [0.003,0.02,0.07,0.16,0.32]:
\n values = self.GetReflectanceValues( Cm, Ch, Bm )
\n CSV_Out = CSV_Out + values
\n groupByHemoglobin.append(values)
\n groupByMelanin.append(groupByHemoglobin)
\n groupByBlend.append(groupByMelanin)
\n CSV_Out_File = open(\"AbsorptionValues.csv\",\"w+\")
\n CSV_Out_File.write(CSV_Out)
\n CSV_Out_File.close()
\n return groupByBlend<\/p>\n
def GetReflectanceValues(self, Cm, Ch, Bm):
\n CSV_Out = \"\"
\n wavelengths = range(380,790,10)
\n for nm in wavelengths:<\/p>\n
# First layer absorption - Epidermis
\n SAV_eumelanin_L = 6.6 * pow(10,10) * pow(nm,-3.33)
\n SAV_pheomelanin_L = 2.9 * pow(10,14) * pow(nm,-4.75)
\n epidermal_hemoglobin_fraction = Ch * 0.25<\/p>\n
# baseline - used in both layers
\n baselineSkinAbsorption_L = 0.0244 + 8.54 * pow(10,-(nm-220)\/123)<\/p>\n
# Epidermis Absorption Coefficient:
\n epidermis = Cm * ((Bm * SAV_eumelanin_L) +((1 - Bm ) * SAV_pheomelanin_L)) + ((1 - epidermal_hemoglobin_fraction) * baselineSkinAbsorption_L)<\/p>\n
# Second layer absorption - Dermis
\n gammaOxy = self.O2Hb[int(nm)]
\n gammaDeoxy = self.Hb[int(nm)]
\n A = 0.75 * gammaOxy + (1 - 0.75) * gammaDeoxy
\n dermis = Ch * (A + (( 1 - Ch) * baselineSkinAbsorption_L))<\/p>\n
# Scattering coefficients
\n scattering_epidermis = pow(14.74 * nm, -0.22) + 2.2 * pow(10,11) * pow(nm, -4)
\n scattering_dermis = scattering_epidermis * 0.5<\/p>\n
reflectance = self.MonteCarlo(epidermis, scattering_epidermis, dermis, scattering_dermis, nm)
\n print(f\"Done on {nm} >> ({Cm}, {Ch}, {Bm}) = {reflectance}\")<\/p>\n
CSV_Out = CSV_Out + f\"{nm},{Cm},{Ch},{Bm},{SAV_eumelanin_L},{SAV_pheomelanin_L},{baselineSkinAbsorption_L},{epidermis},{scattering_epidermis},{scattering_dermis},{dermis},{reflectance}\\n\"
\n return CSV_Out<\/p>\n
def MonteCarlo (self, epi_mua, epi_mus, derm_mua, derm_mus, nm):
\n # These are our Monte Carlo Light Transport Variables that don't change
\n Nbins = 1000
\n Nbinsp1 = 1001
\n PI = 3.1415926
\n LIGHTSPEED = 2.997925 * pow(10,10)
\n ALIVE = 1
\n DEAD = 0
\n THRESHOLD = 0.01
\n CHANCE = 0.1
\n COS90D = 1 * pow(10,-6)
\n ONE_MINUS_COSZERO = 1 * pow(10,-12)
\n COSZERO = 1.0 - 1.0e-12 # cosine of about 1e-6 rad
\n g = 0.9
\n nt = 1.33 # Index of refraction
\n epidermis_thickness = 0.25<\/p>\n
x = 0.0
\n y = 0.0
\n z = 0.0 # photon position<\/p>\n
ux = 0.0
\n uy = 0.0
\n uz = 0.0 # photon trajectory as cosines<\/p>\n
uxx = 0.0
\n uyy = 0.0
\n uzz = 0.0 # temporary values used during SPIN<\/p>\n
s = 0.0 # step sizes. s = -log(RND)\/mus [cm]
\n costheta = 0.0 # cos(theta)
\n sintheta = 0.0 # sin(theta)
\n cospsi = 0.0 # cos(psi)
\n sinpsi = 0.0 # sin(psi)
\n psi = 0.0 # azimuthal angle
\n i_photon = 0.0 # current photon
\n W = 0.0 # photon weight
\n absorb = 0.0 # weighted deposited in a step due to absorption
\n photon_status = 0.0 # flag = ALIVE=1 or DEAD=0
\n ReflBin = [None] * Nbinsp1 #bin to store weights of escaped photos for reflectivity
\n epi_albedo = epi_mus\/(epi_mus + epi_mua) # albedo of tissue
\n derm_albedo = derm_mus\/(derm_mus + derm_mua) # albedo of tissue
\n Nphotons = 1000 # number of photons in simulation
\n NR = Nbins # number of radial positions
\n radial_size = 2.5 # maximum radial size
\n r = 0.0 # radial position
\n dr = radial_size\/NR; # cm, radial bin size
\n ir = 0 # index to radial position
\n shellvolume = 0.0 # volume of shell at radial position r
\n CNT = 0.0 # total count of photon weight summed over all bins
\n rnd = 0.0 # assigned random value 0-1
\n u = 0.0
\n temp = 0.0 # dummy variables<\/p>\n
# Inits
\n random.seed(0)
\n RandomNum = random.random()
\n for i in range(NR+1):
\n ReflBin[i] = 0<\/p>\n
while True:
\n i_photon = i_photon + 1<\/p>\n
W = 1.0
\n photon_status = ALIVE<\/p>\n
x= 0
\n y = 0
\n z = 0<\/p>\n
#Randomly set photon trajectory to yield an isotropic source.
\n costheta = 2.0 * random.random() - 1.0<\/p>\n
sintheta = math.sqrt(1.0 - costheta*costheta)
\n psi = 2.0 * PI * random.random()
\n ux = sintheta * math.cos(psi)
\n uy = sintheta * math.sin(psi)
\n uz = (abs(costheta)) # on the first step we want to head down, into the tissue, so > 0<\/p>\n
# Propagate one photon until it dies as determined by ROULETTE.
\n # or if it reaches the surface again
\n it = 0
\n max_iterations = 100000 # to help avoid infinite loops in case we do something wrong<\/p>\n
# we'll hit epidermis first, so set mua\/mus to those scattering\/absorption values
\n mua = epi_mua
\n mus = epi_mus
\n albedo = epi_albedo
\n while True:
\n it = it + 1
\n rnd = random.random()
\n while rnd <= 0.0: # make sure it is > 0.0
\n rnd = random.random()
\n s = -math.log(rnd)\/(mua + mus)
\n x = x + (s * ux)
\n y = y + (s * uy)
\n z = z + (s * uz)<\/p>\n
if uz < 0:
\n # calculate partial step to reach boundary surface
\n s1 = abs(z\/uz)
\n # move back
\n x = x - (s * ux)
\n y = y - (s * uy)
\n z = z - (s * uz)
\n # take partial step
\n x = x + (s1 * ux)
\n y = y + (s1 * uy)
\n z = z + (s1 * uz)<\/p>\n
# photon is now at the surface boundary, figure out how much escaped and how much was reflected
\n internal_reflectance = self.RFresnel(1.0,nt, -uz )<\/p>\n
#Add weighted reflectance of escpaed photon to reflectance bin
\n external_reflectance = 1 - internal_reflectance
\n r = math.sqrt(x*x + y*y)
\n ir = (r\/dr)
\n if ir >= NR:
\n ir = NR
\n ReflBin[int(ir)] = ReflBin[int(ir)] + (W * external_reflectance)<\/p>\n
# Bounce the photon back into the skin
\n W = internal_reflectance * W
\n uz = -uz
\n x = (s-s1) * ux
\n y = (s-s1) * uy
\n z = (s-s1) * uz<\/p>\n
# check if we have passed into the second layer, or the first
\n if z <= epidermis_thickness:
\n mua = epi_mua
\n mus = epi_mus
\n albedo = epi_albedo
\n else:
\n mua = derm_mua
\n mus = derm_mus
\n albedo = derm_albedo<\/p>\n
''' DROP '''
\n absorb = W*(1 - albedo)
\n W = W - absorb<\/p>\n
''' SPIN '''
\n # Sample for costheta
\n rnd = random.random()
\n if (g == 0.0):
\n costheta = 2.0*rnd - 1.0
\n else:
\n temp = (1.0 - g*g)\/(1.0 - g + 2*g*rnd)
\n costheta = (1.0 + g*g - temp*temp)\/(2.0*g)
\n sintheta = math.sqrt(1.0 - costheta*costheta) <\/p>\n
# Sample psi.
\n psi = 2.0*PI*random.random()
\n cospsi = math.cos(psi)
\n if (psi < PI):
\n sinpsi = math.sqrt(1.0 - cospsi*cospsi)
\n else:
\n sinpsi = -math.sqrt(1.0 - cospsi*cospsi)<\/p>\n
# New trajectory.
\n if (1 - abs(uz) <= ONE_MINUS_COSZERO) : # close to perpendicular.
\n uxx = sintheta * cospsi
\n uyy = sintheta * sinpsi
\n uzz = costheta * self.SIGN(uz) # SIGN() is faster than division.
\n else: # usually use this option
\n temp = math.sqrt(1.0 - uz * uz)
\n uxx = sintheta * (ux * uz * cospsi - uy * sinpsi) \/ temp + ux * costheta
\n uyy = sintheta * (uy * uz * cospsi + ux * sinpsi) \/ temp + uy * costheta
\n uzz = -sintheta * cospsi * temp + uz * costheta<\/p>\n
# Update trajectory
\n ux = uxx
\n uy = uyy
\n uz = uzz<\/p>\n
# Check Roulette
\n if (W < THRESHOLD):
\n if (random.random() <= CHANCE):
\n W = W \/ CHANCE
\n else:
\n photon_status = DEAD
\n if photon_status is DEAD:
\n break
\n if it > max_iterations:
\n break<\/p>\n
if i_photon >= Nphotons:
\n break
\n total_reflection = 0.0
\n for each in range(NR+1):
\n total_reflection = total_reflection + ReflBin[each]\/Nphotons
\n return total_reflection<\/p>\n
def RFresnel(self, n1, n2, cosT1):
\n r = 0.0
\n cosT2 = 0.0
\n COSZERO = 1.0 - 1.0e-12
\n COS90D = 1 * pow(10,-6)
\n if n1 == n2: #matched boundary
\n r = 0.0
\n cosT2 = cosT1
\n elif cosT1 > COSZERO: # normal incident
\n cosT2 = ca1
\n r = (n2-n1)\/(n2+n1)
\n r *= r
\n elif cosT1 < COS90D: # very slant
\n cosT2 = 0.0
\n r = 1.0
\n else: #general
\n sinT1 = math.sqrt(1 - cosT1*cosT1)
\n sinT2 = n1 * sinT1\/n2<\/p>\n
if sinT2 >= 1.0:
\n r = 1.0
\n cosT2 = 0.0
\n else:
\n cosT2 = math.sqrt(1 - sinT2 * sinT2)
\n cosAP = cosT1*cosT2 - sinT1*sinT2
\n cosAM = cosT1*cosT2 + sinT1*sinT2
\n sinAP = sinT1*cosT2 + cosT1*sinT2
\n sinAM = sinT1*cosT2 - cosT1*sinT2
\n r = 0.5 * sinAM * sinAM*(cosAM*cosAM+cosAP*cosAP)\/(sinAP*sinAP*cosAM*cosAM)
\n return r<\/p>\n
def SIGN(self, x):
\n if x >=0:
\n return 1
\n else:
\n return 0 <\/p>\n<\/code>\n\n\n\n
<\/p>\n\n\n\n
skinAbsorption = SkinAbsorption()
\nreflectanceValues = skinAbsorption.Generate()<\/code><\/p>\n\n\n\nimport random
\nimport math<\/code><\/p>\n\n\n\nimport random
\nimport math<\/code>
\n<\/p>\n\n\n\nclass SkinAbsorption:<\/p>\n\n\n\n
O2Hb = {}
\n Hb = {}<\/p>\n\n\n\n
def __init__(self):
\n print(“init SkinAbsorption”)
\n HbFile = open(‘hb.csv’, “r”)
\n O2HbFile = open(‘O2Hb.csv’, “r”)<\/p>\n\n\n\n
HbLines = HbFile.readlines()
\n for line in HbLines:
\n splitLine = line.split(“,”)
\n self.Hb[int(splitLine[0])] = float(splitLine[1].rstrip(“\\n”))<\/p>\n\n\n\n
O2HbLines = O2HbFile.readlines()
\n for line in O2HbLines:
\n splitLine = line.split(“,”)
\n self.O2Hb[int(splitLine[0])] = float(splitLine[1].rstrip(“\\n”))<\/p>\n\n\n\n
HbFile.close()
\n O2HbFile.close()<\/p>\n\n\n\n
def Generate(self):
\n CSV_Out = “Wavelength,Cm, Ch, Bm,Eumelanin,Pheomelanin,Baseline,Epidermis,ScatteringEpi,ScatteringDerm,Dermis,Reflectance\\n”
\n groupByBlend = []
\n for Bm in [0.01,0.5,0.99]:
\n groupByMelanin = []
\n for Cm in [0.002,0.0135,0.0425,0.1,0.185,0.32,0.5]:
\n groupByHemoglobin = []
\n for Ch in [0.003,0.02,0.07,0.16,0.32]:
\n values = self.GetReflectanceValues( Cm, Ch, Bm )
\n CSV_Out = CSV_Out + values
\n groupByHemoglobin.append(values)
\n groupByMelanin.append(groupByHemoglobin)
\n groupByBlend.append(groupByMelanin)
\n CSV_Out_File = open(“AbsorptionValues.csv”,”w+”)
\n CSV_Out_File.write(CSV_Out)
\n CSV_Out_File.close()
\n return groupByBlend<\/p>\n\n\n\n
def GetReflectanceValues(self, Cm, Ch, Bm):
\n CSV_Out = “”
\n wavelengths = range(380,790,10)
\n for nm in wavelengths:<\/p>\n\n\n\n
# First layer absorption – Epidermis
\n SAV_eumelanin_L = 6.6 * pow(10,10) * pow(nm,-3.33)
\n SAV_pheomelanin_L = 2.9 * pow(10,14) * pow(nm,-4.75)
\n epidermal_hemoglobin_fraction = Ch * 0.25<\/p>\n\n\n\n
# baseline – used in both layers
\n baselineSkinAbsorption_L = 0.0244 + 8.54 * pow(10,-(nm-220)\/123)<\/p>\n\n\n\n
# Epidermis Absorption Coefficient:
\n epidermis = Cm * ((Bm * SAV_eumelanin_L) +((1 – Bm ) * SAV_pheomelanin_L)) + ((1 – epidermal_hemoglobin_fraction) * baselineSkinAbsorption_L)<\/p>\n\n\n\n
# Second layer absorption – Dermis
\n gammaOxy = self.O2Hb[int(nm)]
\n gammaDeoxy = self.Hb[int(nm)]
\n A = 0.75 * gammaOxy + (1 – 0.75) * gammaDeoxy
\n dermis = Ch * (A + (( 1 – Ch) * baselineSkinAbsorption_L))<\/p>\n\n\n\n
# Scattering coefficients
\n scattering_epidermis = pow(14.74 * nm, -0.22) + 2.2 * pow(10,11) * pow(nm, -4)
\n scattering_dermis = scattering_epidermis * 0.5<\/p>\n\n\n\n
reflectance = self.MonteCarlo(epidermis, scattering_epidermis, dermis, scattering_dermis, nm)
\n print(f”Done on {nm} >> ({Cm}, {Ch}, {Bm}) = {reflectance}”)<\/p>\n\n\n\n
CSV_Out = CSV_Out + f”{nm},{Cm},{Ch},{Bm},{SAV_eumelanin_L},{SAV_pheomelanin_L},{baselineSkinAbsorption_L},{epidermis},{scattering_epidermis},{scattering_dermis},{dermis},{reflectance}\\n”
\n return CSV_Out<\/p>\n\n\n\n
def MonteCarlo (self, epi_mua, epi_mus, derm_mua, derm_mus, nm):
\n # These are our Monte Carlo Light Transport Variables that don’t change
\n Nbins = 1000
\n Nbinsp1 = 1001
\n PI = 3.1415926
\n LIGHTSPEED = 2.997925 * pow(10,10)
\n ALIVE = 1
\n DEAD = 0
\n THRESHOLD = 0.01
\n CHANCE = 0.1
\n COS90D = 1 * pow(10,-6)
\n ONE_MINUS_COSZERO = 1 * pow(10,-12)
\n COSZERO = 1.0 – 1.0e-12 # cosine of about 1e-6 rad
\n g = 0.9
\n nt = 1.33 # Index of refraction
\n epidermis_thickness = 0.25<\/p>\n\n\n\n
x = 0.0
\n y = 0.0
\n z = 0.0 # photon position<\/p>\n\n\n\n
ux = 0.0
\n uy = 0.0
\n uz = 0.0 # photon trajectory as cosines<\/p>\n\n\n\n
uxx = 0.0
\n uyy = 0.0
\n uzz = 0.0 # temporary values used during SPIN<\/p>\n\n\n\n
s = 0.0 # step sizes. s = -log(RND)\/mus [cm]
\n costheta = 0.0 # cos(theta)
\n sintheta = 0.0 # sin(theta)
\n cospsi = 0.0 # cos(psi)
\n sinpsi = 0.0 # sin(psi)
\n psi = 0.0 # azimuthal angle
\n i_photon = 0.0 # current photon
\n W = 0.0 # photon weight
\n absorb = 0.0 # weighted deposited in a step due to absorption
\n photon_status = 0.0 # flag = ALIVE=1 or DEAD=0
\n ReflBin = [None] * Nbinsp1 #bin to store weights of escaped photos for reflectivity
\n epi_albedo = epi_mus\/(epi_mus + epi_mua) # albedo of tissue
\n derm_albedo = derm_mus\/(derm_mus + derm_mua) # albedo of tissue
\n Nphotons = 1000 # number of photons in simulation
\n NR = Nbins # number of radial positions
\n radial_size = 2.5 # maximum radial size
\n r = 0.0 # radial position
\n dr = radial_size\/NR; # cm, radial bin size
\n ir = 0 # index to radial position
\n shellvolume = 0.0 # volume of shell at radial position r
\n CNT = 0.0 # total count of photon weight summed over all bins
\n rnd = 0.0 # assigned random value 0-1
\n u = 0.0
\n temp = 0.0 # dummy variables<\/p>\n\n\n\n
# Inits
\n random.seed(0)
\n RandomNum = random.random()
\n for i in range(NR+1):
\n ReflBin[i] = 0<\/p>\n\n\n\n
while True:
\n i_photon = i_photon + 1<\/p>\n\n\n\n
W = 1.0
\n photon_status = ALIVE<\/p>\n\n\n\n
x= 0
\n y = 0
\n z = 0<\/p>\n\n\n\n
#Randomly set photon trajectory to yield an isotropic source.
\n costheta = 2.0 * random.random() – 1.0<\/p>\n\n\n\n
sintheta = math.sqrt(1.0 – costheta*costheta)
\n psi = 2.0 * PI * random.random()
\n ux = sintheta * math.cos(psi)
\n uy = sintheta * math.sin(psi)
\n uz = (abs(costheta)) # on the first step we want to head down, into the tissue, so > 0<\/p>\n\n\n\n
# Propagate one photon until it dies as determined by ROULETTE.
\n # or if it reaches the surface again
\n it = 0
\n max_iterations = 100000 # to help avoid infinite loops in case we do something wrong<\/p>\n\n\n\n
# we’ll hit epidermis first, so set mua\/mus to those scattering\/absorption values
\n mua = epi_mua
\n mus = epi_mus
\n albedo = epi_albedo
\n while True:
\n it = it + 1
\n rnd = random.random()
\n while rnd <= 0.0: # make sure it is > 0.0
\n rnd = random.random()
\n s = -math.log(rnd)\/(mua + mus)
\n x = x + (s * ux)
\n y = y + (s * uy)
\n z = z + (s * uz)<\/p>\n\n\n\n
if uz < 0:
\n # calculate partial step to reach boundary surface
\n s1 = abs(z\/uz)
\n # move back
\n x = x – (s * ux)
\n y = y – (s * uy)
\n z = z – (s * uz)
\n # take partial step
\n x = x + (s1 * ux)
\n y = y + (s1 * uy)
\n z = z + (s1 * uz)<\/p>\n\n\n\n
# photon is now at the surface boundary, figure out how much escaped and how much was reflected
\n internal_reflectance = self.RFresnel(1.0,nt, -uz )<\/p>\n\n\n\n
#Add weighted reflectance of escpaed photon to reflectance bin
\n external_reflectance = 1 – internal_reflectance
\n r = math.sqrt(x*x + y*y)
\n ir = (r\/dr)
\n if ir >= NR:
\n ir = NR
\n ReflBin[int(ir)] = ReflBin[int(ir)] + (W * external_reflectance)<\/p>\n\n\n\n
# Bounce the photon back into the skin
\n W = internal_reflectance * W
\n uz = -uz
\n x = (s-s1) * ux
\n y = (s-s1) * uy
\n z = (s-s1) * uz<\/p>\n\n\n\n
# check if we have passed into the second layer, or the first
\n if z <= epidermis_thickness:
\n mua = epi_mua
\n mus = epi_mus
\n albedo = epi_albedo
\n else:
\n mua = derm_mua
\n mus = derm_mus
\n albedo = derm_albedo<\/p>\n\n\n\n
”’ DROP ”’
\n absorb = W*(1 – albedo)
\n W = W – absorb<\/p>\n\n\n\n
”’ SPIN ”’
\n # Sample for costheta
\n rnd = random.random()
\n if (g == 0.0):
\n costheta = 2.0*rnd – 1.0
\n else:
\n temp = (1.0 – g*g)\/(1.0 – g + 2*g*rnd)
\n costheta = (1.0 + g*g – temp*temp)\/(2.0*g)
\n sintheta = math.sqrt(1.0 – costheta*costheta) <\/p>\n\n\n\n
# Sample psi.
\n psi = 2.0*PI*random.random()
\n cospsi = math.cos(psi)
\n if (psi < PI):
\n sinpsi = math.sqrt(1.0 – cospsi*cospsi)
\n else:
\n sinpsi = -math.sqrt(1.0 – cospsi*cospsi)<\/p>\n\n\n\n
# New trajectory.
\n if (1 – abs(uz) <= ONE_MINUS_COSZERO) : # close to perpendicular.
\n uxx = sintheta * cospsi
\n uyy = sintheta * sinpsi
\n uzz = costheta * self.SIGN(uz) # SIGN() is faster than division.
\n else: # usually use this option
\n temp = math.sqrt(1.0 – uz * uz)
\n uxx = sintheta * (ux * uz * cospsi – uy * sinpsi) \/ temp + ux * costheta
\n uyy = sintheta * (uy * uz * cospsi + ux * sinpsi) \/ temp + uy * costheta
\n uzz = -sintheta * cospsi * temp + uz * costheta<\/p>\n\n\n\n
# Update trajectory
\n ux = uxx
\n uy = uyy
\n uz = uzz<\/p>\n\n\n\n
# Check Roulette
\n if (W < THRESHOLD):
\n if (random.random() <= CHANCE):
\n W = W \/ CHANCE
\n else:
\n photon_status = DEAD
\n if photon_status is DEAD:
\n break
\n if it > max_iterations:
\n break<\/p>\n\n\n\n
if i_photon >= Nphotons:
\n break
\n total_reflection = 0.0
\n for each in range(NR+1):
\n total_reflection = total_reflection + ReflBin[each]\/Nphotons
\n return total_reflection<\/p>\n\n\n\n
def RFresnel(self, n1, n2, cosT1):
\n r = 0.0
\n cosT2 = 0.0
\n COSZERO = 1.0 – 1.0e-12
\n COS90D = 1 * pow(10,-6)
\n if n1 == n2: #matched boundary
\n r = 0.0
\n cosT2 = cosT1
\n elif cosT1 > COSZERO: # normal incident
\n cosT2 = ca1
\n r = (n2-n1)\/(n2+n1)
\n r *= r
\n elif cosT1 < COS90D: # very slant
\n cosT2 = 0.0
\n r = 1.0
\n else: #general
\n sinT1 = math.sqrt(1 – cosT1*cosT1)
\n sinT2 = n1 * sinT1\/n2<\/p>\n\n\n\n
![\"\"](\"https:\/\/half4.xyz\/wp-content\/uploads\/2020\/07\/image-36-567x1024.png\")