Part 2 (Practical): The Monte Carlo Simulation

We’ll take what we learned in the previous part and start fleshing out our Monte Carlo (“MC”) simulation.

First, let’s import the random and math classes – we’ll need them.

import random
import math
class SkinAbsorption:
    ....

In our Python script, we’ll add a new function to our class. We’ll call this MonteCarlo. As its inputs, we need these parameters from our spectral work: σ^epidermis_a, σ^epidermis_s, σ^dermis_a, σ^dermis_s and λ

def MonteCarlo (self, epi_mua, epi_mus, derm_mua, derm_mus, nm):

We’ll need two more utility functions. One is the Fresnel function, the second is a quick sign(x) function. Add these to the class, too:

def RFresnel(self, n1, n2, cosT1):
    r = 0.0
    cosT2 = 0.0
    COSZERO = 1.0 - 1.0e-12
    COS90D = 1 * pow(10,-6)
    if n1 == n2: #matched boundary
        r = 0.0
        cosT2 = cosT1
    elif cosT1 > COSZERO:     # normal incident
        cosT2 = ca1
        r = (n2-n1)/(n2+n1)
        r *= r
    elif cosT1 < COS90D:      # very slant
        cosT2 = 0.0
        r = 1.0
    else: #general
        sinT1 = math.sqrt(1 - cosT1*cosT1)
        sinT2 = n1 * sinT1/n2
        
        if sinT2 >= 1.0:
            r = 1.0
            cosT2 = 0.0
        else:
            cosT2 = math.sqrt(1 - sinT2 * sinT2)
            cosAP = cosT1*cosT2 - sinT1*sinT2
            cosAM = cosT1*cosT2 + sinT1*sinT2
            sinAP = sinT1*cosT2 + cosT1*sinT2
            sinAM = sinT1*cosT2 - cosT1*sinT2
            r = 0.5 * sinAM * sinAM*(cosAM*cosAM+cosAP*cosAP)/(sinAP*sinAP*cosAM*cosAM)
    return r

The next step is to add the following values to our Monte Carlo function. These are values we need for our light transport:

# These are our Monte Carlo Light Transport Variables that don't change
Nbins = 1000
Nbinsp1 = 1001
PI = 3.1415926
LIGHTSPEED = exp(2.997925,10)
ALIVE = 1
DEAD = 0
THRESHOLD = 0.01
CHANCE = 0.1
COS90D = exp(1,-6)
ONE_MINUS_COSZERO = exp(1,-12)
COSZERO = 1.0 - 1.0e-12     # cosine of about 1e-6 rad
g = 0.9
nt = 1.4 # Index of refraction
epidermis_thickness = 0.25

Lets now create all the variables we need to change, such as position, trajectory, bins, weight, etc.:

x = 0.0
y = 0.0
z = 0.0 # photon position

ux = 0.0
uy = 0.0
uz = 0.0 # photon trajectory as cosines

uxx = 0.0
uyy = 0.0
uzz = 0.0 # temporary values used during SPIN

s = 0.0 # step sizes. s = -log(RND)/mus [cm] 
costheta = 0.0 # cos(theta) 
sintheta = 0.0 # sin(theta) 
cospsi = 0.0 # cos(psi) 
sinpsi = 0.0 # sin(psi) 
psi = 0.0 # azimuthal angle 
i_photon = 0.0 # current photon
W = 0.0 # photon weight 
absorb = 0.0 # weighted deposited in a step due to absorption 
photon_status = 0.0 # flag = ALIVE=1 or DEAD=0 
ReflBin = [None] * Nbinsp1 #bin to store weights of escaped photos for reflectivity
epi_albedo = epi_mus/(epi_mus + epi_mua) # albedo of tissue
derm_albedo = derm_mus/(derm_mus + derm_mua) # albedo of tissue
Nphotons = 1000 # number of photons in simulation 
NR = Nbins # number of radial positions 
radial_size = 2.5 # maximum radial size 
r = 0.0 # radial position 
dr = radial_size/NR; # cm, radial bin size 
ir = 0 # index to radial position 
shellvolume = 0.0 # volume of shell at radial position r 
CNT = 0.0 # total count of photon weight summed over all bins 
rnd = 0.0 # assigned random value 0-1 
u = 0.0
temp = 0.0 # dummy variables

And the last thing before we begin looping; some initialization:

# Inits
random.seed(0)
RandomNum = random.random()
for i in range(NR+1):
    ReflBin[i] = 0

Now we can begin our iteration code. We’ll start with a simple while loop for number of photons:

while True:
    i_photon = i_photon + 1
    if i_photon >= Nphotons:
        break

Within this loop, we need to do some further initialization for this loop:

W = 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

W = 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

W = 1.0
photon_status = ALIVE

W = 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

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.

Let’s make a new loop within the existing one. The first thing we have to do is add a break condition – for now, we’ll put in a max iteration, but we’ll also add in our terminate functionality.

while True:            
    it = it + 1
    # Check Roulette
    if (W < THRESHOLD):
        if (random.random() <= CHANCE):
            W = W / CHANCE
        else:
            photon_status = DEAD
    if photon_status is DEAD:
        break
    if it > max_iterations:
        break

while True:
it = it + 1
# Check Roulette
if (W < THRESHOLD):
if (random.random() <= CHANCE):
W = W / CHANCE
else:
photon_status = DEAD
if photon_status is DEAD:
break
if it > max_iterations:
break

while True:
it = it + 1
# Check Roulette
if (W < THRESHOLD):
if (random.random() <= CHANCE):
W = W / CHANCE
else:
photon_status = DEAD
if photon_status is DEAD:
break
if it > max_iterations:
break

while True:
it = it + 1
# Check Roulette
if (W < THRESHOLD):
if (random.random() <= CHANCE):
W = W / CHANCE
else:
photon_status = DEAD
if photon_status is DEAD:
break
if it > max_iterations:
break

Above the ‘# Check Roulette‘ line, we’ll flesh out our hop code:

# Check Roulette

# Check Roulette

# Check Roulette

rnd = random.random()
while rnd <= 0.0: # make sure it is > 0.0
    rnd = random.random()
s = -math.log(rnd)/(mua + mus)
x = x + (s * ux)
y = y + (s * uy)
z = z + (s * uz)

rnd = random.random()
while rnd <= 0.0: # make sure it is > 0.0
rnd = random.random()
s = -math.log(rnd)/(mua + mus)
x = x + (s * ux)
y = y + (s * uy)
z = z + (s * uz)

rnd = random.random()
while rnd <= 0.0: # make sure it is > 0.0
rnd = random.random()
s = -math.log(rnd)/(mua + mus)
x = x + (s * ux)
y = y + (s * uy)
z = z + (s * uz)

rnd = random.random()
while rnd <= 0.0: # make sure it is > 0.0
rnd = random.random()
s = -math.log(rnd)/(mua + mus)
x = x + (s * ux)
y = y + (s * uy)
z = z + (s * uz)

On our first ‘hop’, the code is the ‘launch’ 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:

if uz < 0:
    # calculate partial step to reach boundary surface
    s1 = abs(z/uz)
    # move back
    x = x - (s * ux)
    y = y - (s * uy)
    z = z - (s * uz)
    # take partial step
    x = x + (s1 * ux)
    y = y + (s1 * uy)
    z = z + (s1 * uz)

if uz < 0:
# calculate partial step to reach boundary surface
s1 = abs(z/uz)
# move back
x = x - (s * ux)
y = y - (s * uy)
z = z - (s * uz)
# take partial step
x = x + (s1 * ux)
y = y + (s1 * uy)
z = z + (s1 * uz)

if uz < 0:
# calculate partial step to reach boundary surface
s1 = abs(z/uz)
# move back
x = x - (s * ux)
y = y - (s * uy)
z = z - (s * uz)
# take partial step
x = x + (s1 * ux)
y = y + (s1 * uy)
z = z + (s1 * uz)

if uz < 0:
# calculate partial step to reach boundary surface
s1 = abs(z/uz)
# move back
x = x - (s * ux)
y = y - (s * uy)
z = z - (s * uz)
# take partial step
x = x + (s1 * ux)
y = y + (s1 * uy)
z = z + (s1 * uz)

To 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:

# photon is now at the surface boundary, figure out how much escaped and how much was reflected
internal_reflectance = self.RFresnel(1.0,nt, -uz )

#Add weighted reflectance of escpaed photon to reflectance bin
external_reflectance = 1 - internal_reflectance
r = math.sqrt(x*x + y*y)
ir = (r/dr)
if ir >= NR:
   ir = NR  
ReflBin[int(ir)] = ReflBin[int(ir)] + (W * external_reflectance)
						
# Bounce the photon back into the skin
W = internal_reflectance * W
uz = -uz
x = (s-s1) * ux
y = (s-s1) * uy
z = (s-s1) * uz

# photon is now at the surface boundary, figure out how much escaped and how much was reflected
internal_reflectance = self.RFresnel(1.0,nt, -uz )

#Add weighted reflectance of escpaed photon to reflectance bin
external_reflectance = 1 - internal_reflectance
r = math.sqrt(x*x + y*y)
ir = (r/dr)
if ir >= NR:
ir = NR
ReflBin[int(ir)] = ReflBin[int(ir)] + (W * external_reflectance)

# Bounce the photon back into the skin
W = internal_reflectance * W
uz = -uz
x = (s-s1) * ux
y = (s-s1) * uy
z = (s-s1) * uz

# photon is now at the surface boundary, figure out how much escaped and how much was reflected
internal_reflectance = self.RFresnel(1.0,nt, -uz )

# photon is now at the surface boundary, figure out how much escaped and how much was reflected
internal_reflectance = self.RFresnel(1.0,nt, -uz )

#Add weighted reflectance of escpaed photon to reflectance bin
external_reflectance = 1 – internal_reflectance
r = math.sqrt(x*x + y*y)
ir = (r/dr)
if ir >= NR:
ir = NR
ReflBin[int(ir)] = ReflBin[int(ir)] + (W * external_reflectance)

# Bounce the photon back into the skin
W = internal_reflectance * W
uz = -uz
x = (s-s1) * ux
y = (s-s1) * uy
z = (s-s1) * uz

Close of the if-statement. Rather than do a full boundary check between the skin layers, we’re 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.

# check if we have passed into the second layer, or the first
if z <= epidermis_thickness:
   mua = epi_mua
   mus = epi_mus
   albedo = epi_albedo
else:
   mua = derm_mua
   mus = derm_mus
   albedo = derm_albedo

# check if we have passed into the second layer, or the first
if z <= epidermis_thickness:
mua = epi_mua
mus = epi_mus
albedo = epi_albedo
else:
mua = derm_mua
mus = derm_mus
albedo = derm_albedo

# check if we have passed into the second layer, or the first
if z <= epidermis_thickness:
mua = epi_mua
mus = epi_mus
albedo = epi_albedo
else:
mua = derm_mua
mus = derm_mus
albedo = derm_albedo

# check if we have passed into the second layer, or the first
if z <= epidermis_thickness:
mua = epi_mua
mus = epi_mus
albedo = epi_albedo
else:
mua = derm_mua
mus = derm_mus
albedo = derm_albedo

Drop stage now – fairly simple. We just reduce the photon weight by the albedo.

absorb = W*(1 - albedo)
W = W - absorb

absorb = W*(1 - albedo)
W = W - absorb

absorb = W*(1 - albedo)
W = W - absorb

absorb = W*(1 - albedo)
W = W - absorb

Straight into the Spin… :

''' SPIN 
   Scatter photon into new trajectory defined by theta and psi.
   Theta is specified by cos(theta), which is determined 
   based on the Henyey-Greenstein scattering function.
   Convert theta and psi into cosines ux, uy, uz. 
'''
# Sample for costheta 
rnd = random.random()
if (g == 0.0):
   costheta = 2.0*rnd - 1.0
else:
   temp = (1.0 - g*g)/(1.0 - g + 2*g*rnd)
   costheta = (1.0 + g*g - temp*temp)/(2.0*g)
sintheta = math.sqrt(1.0 - costheta*costheta) 

# Sample psi. 
psi = 2.0*PI*random.random()
cospsi = math.cos(psi)
if (psi < PI):
   sinpsi = math.sqrt(1.0 - cospsi*cospsi)
else:
   sinpsi = -math.sqrt(1.0 - cospsi*cospsi)

# New trajectory. 
if (1 - abs(uz) <= ONE_MINUS_COSZERO) :      # close to perpendicular. 
   uxx = sintheta * cospsi
   uyy = sintheta * sinpsi
   uzz = costheta * self.SIGN(uz)   # SIGN() is faster than division. 
else: # usually use this option 
   temp = math.sqrt(1.0 - uz * uz)
   uxx = sintheta * (ux * uz * cospsi - uy * sinpsi) / temp + ux * costheta
   uyy = sintheta * (uy * uz * cospsi + ux * sinpsi) / temp + uy * costheta
   uzz = -sintheta * cospsi * temp + uz * costheta

# Update trajectory 
ux = uxx
uy = uyy
uz = uzz

''' SPIN
Scatter photon into new trajectory defined by theta and psi.
Theta is specified by cos(theta), which is determined
based on the Henyey-Greenstein scattering function.
Convert theta and psi into cosines ux, uy, uz.
'''
# Sample for costheta
rnd = random.random()
if (g == 0.0):
costheta = 2.0*rnd - 1.0
else:
temp = (1.0 - g*g)/(1.0 - g + 2*g*rnd)
costheta = (1.0 + g*g - temp*temp)/(2.0*g)
sintheta = math.sqrt(1.0 - costheta*costheta)

# Sample psi.
psi = 2.0*PI*random.random()
cospsi = math.cos(psi)
if (psi < PI):
sinpsi = math.sqrt(1.0 - cospsi*cospsi)
else:
sinpsi = -math.sqrt(1.0 - cospsi*cospsi)

# New trajectory.
if (1 - abs(uz) <= ONE_MINUS_COSZERO) : # close to perpendicular.
uxx = sintheta * cospsi
uyy = sintheta * sinpsi
uzz = costheta * self.SIGN(uz) # SIGN() is faster than division.
else: # usually use this option
temp = math.sqrt(1.0 - uz * uz)
uxx = sintheta * (ux * uz * cospsi - uy * sinpsi) / temp + ux * costheta
uyy = sintheta * (uy * uz * cospsi + ux * sinpsi) / temp + uy * costheta
uzz = -sintheta * cospsi * temp + uz * costheta

# Update trajectory
ux = uxx
uy = uyy
uz = uzz

''' SPIN
Scatter photon into new trajectory defined by theta and psi.
Theta is specified by cos(theta), which is determined
based on the Henyey-Greenstein scattering function.
Convert theta and psi into cosines ux, uy, uz.
'''
# Sample for costheta
rnd = random.random()
if (g == 0.0):
costheta = 2.0*rnd - 1.0
else:
temp = (1.0 - g*g)/(1.0 - g + 2*g*rnd)
costheta = (1.0 + g*g - temp*temp)/(2.0*g)
sintheta = math.sqrt(1.0 - costheta*costheta)

''' SPIN
Scatter photon into new trajectory defined by theta and psi.
Theta is specified by cos(theta), which is determined
based on the Henyey-Greenstein scattering function.
Convert theta and psi into cosines ux, uy, uz.
'''
# Sample for costheta
rnd = random.random()
if (g == 0.0):
costheta = 2.0*rnd - 1.0
else:
temp = (1.0 - g*g)/(1.0 - g + 2*g*rnd)
costheta = (1.0 + g*g - temp*temp)/(2.0*g)
sintheta = math.sqrt(1.0 - costheta*costheta)

# Sample psi.
psi = 2.0*PI*random.random()
cospsi = math.cos(psi)
if (psi < PI):
sinpsi = math.sqrt(1.0 – cospsi*cospsi)
else:
sinpsi = -math.sqrt(1.0 – cospsi*cospsi)

# New trajectory.
if (1 – abs(uz) <= ONE_MINUS_COSZERO) : # close to perpendicular.
uxx = sintheta * cospsi
uyy = sintheta * sinpsi
uzz = costheta * self.SIGN(uz) # SIGN() is faster than division.
else: # usually use this option
temp = math.sqrt(1.0 – uz * uz)
uxx = sintheta * (ux * uz * cospsi – uy * sinpsi) / temp + ux * costheta
uyy = sintheta * (uy * uz * cospsi + ux * sinpsi) / temp + uy * costheta
uzz = -sintheta * cospsi * temp + uz * costheta

# Update trajectory
ux = uxx
uy = uyy
uz = uzz

The terminate code we wrote towards the beginning should sit below all of the above, being the final stage.

Calculating Reflectance

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:

total_reflection = 0.0
   for each in range(NR+1):
      total_reflection = total_reflection + ReflBin[each]/Nphotons
   return total_reflection

total_reflection = 0.0
for each in range(NR+1):
total_reflection = total_reflection + ReflBin[each]/Nphotons
return total_reflection

total_reflection = 0.0
for each in range(NR+1):
total_reflection = total_reflection + ReflBin[each]/Nphotons
return total_reflection

total_reflection = 0.0
for each in range(NR+1):
total_reflection = total_reflection + ReflBin[each]/Nphotons
return total_reflection

Now we need up update our GetReflectanceValues() function so that it can manage the data returned. For debugging purposes, we’ll put the reflectance values on the end of our CSV debug code, like so:

GetReflectanceValues

GetReflectanceValues

[....]
# Scattering coefficients
scattering_epidermis = pow(14.74 * nm, -0.22) + 2.2 * pow(10,11) * pow(nm, -4)
scattering_dermis = scattering_epidermis * 0.5
            
reflectance = self.MonteCarlo(epidermis, scattering_epidermis, dermis, scattering_dermis, nm)
            
CSV_Out = CSV_Out + f"{nm},{SAV_eumelanin_L},{SAV_pheomelanin_L},{baselineSkinAbsorption_L},{epidermis},{scattering_epidermis},{scattering_dermis},{dermis},{reflectance}\n"
 return CSV_Out
[...]

[....]
# Scattering coefficients
scattering_epidermis = pow(14.74 * nm, -0.22) + 2.2 * pow(10,11) * pow(nm, -4)
scattering_dermis = scattering_epidermis * 0.5

reflectance = self.MonteCarlo(epidermis, scattering_epidermis, dermis, scattering_dermis, nm)

CSV_Out = CSV_Out + f"{nm},{SAV_eumelanin_L},{SAV_pheomelanin_L},{baselineSkinAbsorption_L},{epidermis},{scattering_epidermis},{scattering_dermis},{dermis},{reflectance}\n"
return CSV_Out
[...]

[....]
# Scattering coefficients
scattering_epidermis = pow(14.74 * nm, -0.22) + 2.2 * pow(10,11) * pow(nm, -4)
scattering_dermis = scattering_epidermis * 0.5

[....]
# Scattering coefficients
scattering_epidermis = pow(14.74 * nm, -0.22) + 2.2 * pow(10,11) * pow(nm, -4)
scattering_dermis = scattering_epidermis * 0.5

reflectance = self.MonteCarlo(epidermis, scattering_epidermis, dermis, scattering_dermis, nm)

CSV_Out = CSV_Out + f”{nm},{SAV_eumelanin_L},{SAV_pheomelanin_L},{baselineSkinAbsorption_L},{epidermis},{scattering_epidermis},{scattering_dermis},{dermis},{reflectance}\n”
return CSV_Out
[…]

Also don’t forget to update the headers in Generate()…

Generate()

Generate()

Generate()

def Generate(self):
    CSV_Out = "Wavelength,Eumelanin,Pheomelanin,Baseline,Epidermis,ScatteringEpi,ScatteringDerm,Dermis,Reflectance\n"

def Generate(self):
CSV_Out = "Wavelength,Eumelanin,Pheomelanin,Baseline,Epidermis,ScatteringEpi,ScatteringDerm,Dermis,Reflectance\n"


def Generate(self):
CSV_Out = "Wavelength,Eumelanin,Pheomelanin,Baseline,Epidermis,ScatteringEpi,ScatteringDerm,Dermis,Reflectance\n"


def Generate(self):
CSV_Out = "Wavelength,Eumelanin,Pheomelanin,Baseline,Epidermis,ScatteringEpi,ScatteringDerm,Dermis,Reflectance\n"


Debugging

Give the code a run and import your CSV into Excel again.This may take some time, it takes a good 15 mins on a Core i9-9900k. I haven’t implemented multi-threading for this tutorial as it complicates what is already a heavy topic. If you’re confident modifying this to be multi-threaded, feel free to do so.

This may take some time
This may take some time
This may take some time

This may take some time
This may take some time
This may take some time

This may take some time
This may take some time
This may take some time

This may take some time
This may take some time

This may take some time

The final code:

import random
import math

class SkinAbsorption:

    O2Hb = {}
    Hb = {}
    
    def __init__(self):
        print("init SkinAbsorption")
        HbFile = open('hb.csv', "r")
        O2HbFile = open('O2Hb.csv', "r")

        HbLines = HbFile.readlines()
        for line in HbLines:
            splitLine = line.split(",")
            self.Hb[int(splitLine[0])] = float(splitLine[1].rstrip("\n"))

        O2HbLines = O2HbFile.readlines()
        for line in O2HbLines:
            splitLine = line.split(",")
            self.O2Hb[int(splitLine[0])] = float(splitLine[1].rstrip("\n"))
            
        HbFile.close()
        O2HbFile.close()

    def Generate(self):
        CSV_Out = "Wavelength,Cm, Ch, Bm,Eumelanin,Pheomelanin,Baseline,Epidermis,ScatteringEpi,ScatteringDerm,Dermis,Reflectance\n"
        groupByBlend = []
        for Bm in [0.01,0.5,0.99]:
            groupByMelanin = []
            for Cm in [0.002,0.0135,0.0425,0.1,0.185,0.32,0.5]:
                groupByHemoglobin = []
                for Ch in [0.003,0.02,0.07,0.16,0.32]:
                    values = self.GetReflectanceValues( Cm, Ch, Bm )
                    CSV_Out = CSV_Out + values
                    groupByHemoglobin.append(values)
                groupByMelanin.append(groupByHemoglobin)
            groupByBlend.append(groupByMelanin)
        CSV_Out_File = open("AbsorptionValues.csv","w+")
        CSV_Out_File.write(CSV_Out)
        CSV_Out_File.close() 
        return groupByBlend
    
    def GetReflectanceValues(self, Cm, Ch, Bm):
        CSV_Out = ""
        wavelengths = range(380,790,10)
        for nm in wavelengths:
            
            # First layer absorption - Epidermis
            SAV_eumelanin_L = 6.6 * pow(10,10) * pow(nm,-3.33)
            SAV_pheomelanin_L = 2.9 * pow(10,14) * pow(nm,-4.75)
            epidermal_hemoglobin_fraction = Ch * 0.25
            
            # baseline - used in both layers
            baselineSkinAbsorption_L = 0.0244 + 8.54 * pow(10,-(nm-220)/123)

            # Epidermis Absorption Coefficient:          
            epidermis = Cm * ((Bm * SAV_eumelanin_L) +((1 -  Bm ) * SAV_pheomelanin_L)) + ((1 - epidermal_hemoglobin_fraction) * baselineSkinAbsorption_L)
            
            # Second layer absorption - Dermis
            gammaOxy = self.O2Hb[int(nm)]
            gammaDeoxy = self.Hb[int(nm)]
            A = 0.75 * gammaOxy + (1 - 0.75) * gammaDeoxy      
            dermis = Ch * (A + (( 1 - Ch) * baselineSkinAbsorption_L))
            
            # Scattering coefficients
            scattering_epidermis = pow(14.74 * nm, -0.22) + 2.2 * pow(10,11) * pow(nm, -4)
            scattering_dermis = scattering_epidermis * 0.5
            
            reflectance = self.MonteCarlo(epidermis, scattering_epidermis, dermis, scattering_dermis, nm)
            print(f"Done on {nm} >> ({Cm}, {Ch}, {Bm}) = {reflectance}")
            
            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"
        return CSV_Out

    def MonteCarlo (self, epi_mua, epi_mus, derm_mua, derm_mus, nm):
        # These are our Monte Carlo Light Transport Variables that don't change
        Nbins = 1000
        Nbinsp1 = 1001
        PI = 3.1415926
        LIGHTSPEED = 2.997925 * pow(10,10)
        ALIVE = 1
        DEAD = 0
        THRESHOLD = 0.01
        CHANCE = 0.1
        COS90D = 1 * pow(10,-6)
        ONE_MINUS_COSZERO = 1 * pow(10,-12)
        COSZERO = 1.0 - 1.0e-12     # cosine of about 1e-6 rad
        g = 0.9
        nt = 1.33 # Index of refraction
        epidermis_thickness = 0.25

        x = 0.0
        y = 0.0
        z = 0.0 # photon position

        ux = 0.0
        uy = 0.0
        uz = 0.0 # photon trajectory as cosines

        uxx = 0.0
        uyy = 0.0
        uzz = 0.0 # temporary values used during SPIN

        s = 0.0 # step sizes. s = -log(RND)/mus [cm] 
        costheta = 0.0 # cos(theta) 
        sintheta = 0.0 # sin(theta) 
        cospsi = 0.0 # cos(psi) 
        sinpsi = 0.0 # sin(psi) 
        psi = 0.0 # azimuthal angle 
        i_photon = 0.0 # current photon
        W = 0.0 # photon weight 
        absorb = 0.0 # weighted deposited in a step due to absorption 
        photon_status = 0.0 # flag = ALIVE=1 or DEAD=0 
        ReflBin = [None] * Nbinsp1 #bin to store weights of escaped photos for reflectivity
        epi_albedo = epi_mus/(epi_mus + epi_mua) # albedo of tissue
        derm_albedo = derm_mus/(derm_mus + derm_mua) # albedo of tissue
        Nphotons = 1000 # number of photons in simulation 
        NR = Nbins # number of radial positions 
        radial_size = 2.5 # maximum radial size 
        r = 0.0 # radial position 
        dr = radial_size/NR; # cm, radial bin size 
        ir = 0 # index to radial position 
        shellvolume = 0.0 # volume of shell at radial position r 
        CNT = 0.0 # total count of photon weight summed over all bins 
        rnd = 0.0 # assigned random value 0-1 
        u = 0.0
        temp = 0.0 # dummy variables

        # Inits
        random.seed(0)
        RandomNum = random.random()
        for i in range(NR+1):
            ReflBin[i] = 0

        while True:
            i_photon = i_photon + 1

            W = 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
            while True:            
                it = it + 1
                rnd = random.random()
                while rnd <= 0.0: # make sure it is > 0.0
                    rnd = random.random()
                s = -math.log(rnd)/(mua + mus)
                x = x + (s * ux)
                y = y + (s * uy)
                z = z + (s * uz)

                if uz < 0:
                    # calculate partial step to reach boundary surface
                    s1 = abs(z/uz)
                    # move back
                    x = x - (s * ux)
                    y = y - (s * uy)
                    z = z - (s * uz)
                    # take partial step
                    x = x + (s1 * ux)
                    y = y + (s1 * uy)
                    z = z + (s1 * uz)

                    # photon is now at the surface boundary, figure out how much escaped and how much was reflected
                    internal_reflectance = self.RFresnel(1.0,nt, -uz )

                    #Add weighted reflectance of escpaed photon to reflectance bin
                    external_reflectance = 1 - internal_reflectance
                    r = math.sqrt(x*x + y*y)
                    ir = (r/dr)
                    if ir >= NR:
                       ir = NR  
                    ReflBin[int(ir)] = ReflBin[int(ir)] + (W * external_reflectance)
                                                                    
                    # Bounce the photon back into the skin
                    W = internal_reflectance * W
                    uz = -uz
                    x = (s-s1) * ux
                    y = (s-s1) * uy
                    z = (s-s1) * uz

                # check if we have passed into the second layer, or the first
                if z <= epidermis_thickness:
                   mua = epi_mua
                   mus = epi_mus
                   albedo = epi_albedo
                else:
                   mua = derm_mua
                   mus = derm_mus
                   albedo = derm_albedo

                ''' DROP '''
                absorb = W*(1 - albedo)
                W = W - absorb

                ''' SPIN '''
                # Sample for costheta 
                rnd = random.random()
                if (g == 0.0):
                   costheta = 2.0*rnd - 1.0
                else:
                   temp = (1.0 - g*g)/(1.0 - g + 2*g*rnd)
                   costheta = (1.0 + g*g - temp*temp)/(2.0*g)
                sintheta = math.sqrt(1.0 - costheta*costheta) 

                # Sample psi. 
                psi = 2.0*PI*random.random()
                cospsi = math.cos(psi)
                if (psi < PI):
                   sinpsi = math.sqrt(1.0 - cospsi*cospsi)
                else:
                   sinpsi = -math.sqrt(1.0 - cospsi*cospsi)

                # New trajectory. 
                if (1 - abs(uz) <= ONE_MINUS_COSZERO) :      # close to perpendicular. 
                   uxx = sintheta * cospsi
                   uyy = sintheta * sinpsi
                   uzz = costheta * self.SIGN(uz)   # SIGN() is faster than division. 
                else: # usually use this option 
                   temp = math.sqrt(1.0 - uz * uz)
                   uxx = sintheta * (ux * uz * cospsi - uy * sinpsi) / temp + ux * costheta
                   uyy = sintheta * (uy * uz * cospsi + ux * sinpsi) / temp + uy * costheta
                   uzz = -sintheta * cospsi * temp + uz * costheta

                # Update trajectory 
                ux = uxx
                uy = uyy
                uz = uzz

            
                # Check Roulette
                if (W < THRESHOLD):
                    if (random.random() <= CHANCE):
                        W = W / CHANCE
                    else:
                        photon_status = DEAD
                if photon_status is DEAD:
                    break
                if it > max_iterations:
                    break
            
            if i_photon >= Nphotons:
                break
        total_reflection = 0.0
        for each in range(NR+1):
            total_reflection = total_reflection + ReflBin[each]/Nphotons
        return total_reflection

    def RFresnel(self, n1, n2, cosT1):
        r = 0.0
        cosT2 = 0.0
        COSZERO = 1.0 - 1.0e-12
        COS90D = 1 * pow(10,-6)
        if n1 == n2: #matched boundary
            r = 0.0
            cosT2 = cosT1
        elif cosT1 > COSZERO:     # normal incident
            cosT2 = ca1
            r = (n2-n1)/(n2+n1)
            r *= r
        elif cosT1 < COS90D:      # very slant
            cosT2 = 0.0
            r = 1.0
        else: #general
            sinT1 = math.sqrt(1 - cosT1*cosT1)
            sinT2 = n1 * sinT1/n2
            
            if sinT2 >= 1.0:
                r = 1.0
                cosT2 = 0.0
            else:
                cosT2 = math.sqrt(1 - sinT2 * sinT2)
                cosAP = cosT1*cosT2 - sinT1*sinT2
                cosAM = cosT1*cosT2 + sinT1*sinT2
                sinAP = sinT1*cosT2 + cosT1*sinT2
                sinAM = sinT1*cosT2 - cosT1*sinT2
                r = 0.5 * sinAM * sinAM*(cosAM*cosAM+cosAP*cosAP)/(sinAP*sinAP*cosAM*cosAM)
        return r

    def SIGN(self, x):
        if x >=0:
            return 1
        else:
            return 0 

skinAbsorption = SkinAbsorption()
reflectanceValues = skinAbsorption.Generate()

import random
import math

class SkinAbsorption:

O2Hb = {}
Hb = {}

def __init__(self):
print("init SkinAbsorption")
HbFile = open('hb.csv', "r")
O2HbFile = open('O2Hb.csv', "r")

HbLines = HbFile.readlines()
for line in HbLines:
splitLine = line.split(",")
self.Hb[int(splitLine[0])] = float(splitLine[1].rstrip("\n"))

O2HbLines = O2HbFile.readlines()
for line in O2HbLines:
splitLine = line.split(",")
self.O2Hb[int(splitLine[0])] = float(splitLine[1].rstrip("\n"))

HbFile.close()
O2HbFile.close()

def Generate(self):
CSV_Out = "Wavelength,Cm, Ch, Bm,Eumelanin,Pheomelanin,Baseline,Epidermis,ScatteringEpi,ScatteringDerm,Dermis,Reflectance\n"
groupByBlend = []
for Bm in [0.01,0.5,0.99]:
groupByMelanin = []
for Cm in [0.002,0.0135,0.0425,0.1,0.185,0.32,0.5]:
groupByHemoglobin = []
for Ch in [0.003,0.02,0.07,0.16,0.32]:
values = self.GetReflectanceValues( Cm, Ch, Bm )
CSV_Out = CSV_Out + values
groupByHemoglobin.append(values)
groupByMelanin.append(groupByHemoglobin)
groupByBlend.append(groupByMelanin)
CSV_Out_File = open("AbsorptionValues.csv","w+")
CSV_Out_File.write(CSV_Out)
CSV_Out_File.close()
return groupByBlend

def GetReflectanceValues(self, Cm, Ch, Bm):
CSV_Out = ""
wavelengths = range(380,790,10)
for nm in wavelengths:

# First layer absorption - Epidermis
SAV_eumelanin_L = 6.6 * pow(10,10) * pow(nm,-3.33)
SAV_pheomelanin_L = 2.9 * pow(10,14) * pow(nm,-4.75)
epidermal_hemoglobin_fraction = Ch * 0.25

# baseline - used in both layers
baselineSkinAbsorption_L = 0.0244 + 8.54 * pow(10,-(nm-220)/123)

# Epidermis Absorption Coefficient:
epidermis = Cm * ((Bm * SAV_eumelanin_L) +((1 - Bm ) * SAV_pheomelanin_L)) + ((1 - epidermal_hemoglobin_fraction) * baselineSkinAbsorption_L)

# Second layer absorption - Dermis
gammaOxy = self.O2Hb[int(nm)]
gammaDeoxy = self.Hb[int(nm)]
A = 0.75 * gammaOxy + (1 - 0.75) * gammaDeoxy
dermis = Ch * (A + (( 1 - Ch) * baselineSkinAbsorption_L))

# Scattering coefficients
scattering_epidermis = pow(14.74 * nm, -0.22) + 2.2 * pow(10,11) * pow(nm, -4)
scattering_dermis = scattering_epidermis * 0.5

reflectance = self.MonteCarlo(epidermis, scattering_epidermis, dermis, scattering_dermis, nm)
print(f"Done on {nm} >> ({Cm}, {Ch}, {Bm}) = {reflectance}")

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"
return CSV_Out

def MonteCarlo (self, epi_mua, epi_mus, derm_mua, derm_mus, nm):
# These are our Monte Carlo Light Transport Variables that don't change
Nbins = 1000
Nbinsp1 = 1001
PI = 3.1415926
LIGHTSPEED = 2.997925 * pow(10,10)
ALIVE = 1
DEAD = 0
THRESHOLD = 0.01
CHANCE = 0.1
COS90D = 1 * pow(10,-6)
ONE_MINUS_COSZERO = 1 * pow(10,-12)
COSZERO = 1.0 - 1.0e-12 # cosine of about 1e-6 rad
g = 0.9
nt = 1.33 # Index of refraction
epidermis_thickness = 0.25

x = 0.0
y = 0.0
z = 0.0 # photon position

ux = 0.0
uy = 0.0
uz = 0.0 # photon trajectory as cosines

uxx = 0.0
uyy = 0.0
uzz = 0.0 # temporary values used during SPIN

s = 0.0 # step sizes. s = -log(RND)/mus [cm]
costheta = 0.0 # cos(theta)
sintheta = 0.0 # sin(theta)
cospsi = 0.0 # cos(psi)
sinpsi = 0.0 # sin(psi)
psi = 0.0 # azimuthal angle
i_photon = 0.0 # current photon
W = 0.0 # photon weight
absorb = 0.0 # weighted deposited in a step due to absorption
photon_status = 0.0 # flag = ALIVE=1 or DEAD=0
ReflBin = [None] * Nbinsp1 #bin to store weights of escaped photos for reflectivity
epi_albedo = epi_mus/(epi_mus + epi_mua) # albedo of tissue
derm_albedo = derm_mus/(derm_mus + derm_mua) # albedo of tissue
Nphotons = 1000 # number of photons in simulation
NR = Nbins # number of radial positions
radial_size = 2.5 # maximum radial size
r = 0.0 # radial position
dr = radial_size/NR; # cm, radial bin size
ir = 0 # index to radial position
shellvolume = 0.0 # volume of shell at radial position r
CNT = 0.0 # total count of photon weight summed over all bins
rnd = 0.0 # assigned random value 0-1
u = 0.0
temp = 0.0 # dummy variables

# Inits
random.seed(0)
RandomNum = random.random()
for i in range(NR+1):
ReflBin[i] = 0

while True:
i_photon = i_photon + 1

W = 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
while True:
it = it + 1
rnd = random.random()
while rnd <= 0.0: # make sure it is > 0.0
rnd = random.random()
s = -math.log(rnd)/(mua + mus)
x = x + (s * ux)
y = y + (s * uy)
z = z + (s * uz)

if uz < 0:
# calculate partial step to reach boundary surface
s1 = abs(z/uz)
# move back
x = x - (s * ux)
y = y - (s * uy)
z = z - (s * uz)
# take partial step
x = x + (s1 * ux)
y = y + (s1 * uy)
z = z + (s1 * uz)

# photon is now at the surface boundary, figure out how much escaped and how much was reflected
internal_reflectance = self.RFresnel(1.0,nt, -uz )

#Add weighted reflectance of escpaed photon to reflectance bin
external_reflectance = 1 - internal_reflectance
r = math.sqrt(x*x + y*y)
ir = (r/dr)
if ir >= NR:
ir = NR
ReflBin[int(ir)] = ReflBin[int(ir)] + (W * external_reflectance)

# Bounce the photon back into the skin
W = internal_reflectance * W
uz = -uz
x = (s-s1) * ux
y = (s-s1) * uy
z = (s-s1) * uz

# check if we have passed into the second layer, or the first
if z <= epidermis_thickness:
mua = epi_mua
mus = epi_mus
albedo = epi_albedo
else:
mua = derm_mua
mus = derm_mus
albedo = derm_albedo

''' DROP '''
absorb = W*(1 - albedo)
W = W - absorb

''' SPIN '''
# Sample for costheta
rnd = random.random()
if (g == 0.0):
costheta = 2.0*rnd - 1.0
else:
temp = (1.0 - g*g)/(1.0 - g + 2*g*rnd)
costheta = (1.0 + g*g - temp*temp)/(2.0*g)
sintheta = math.sqrt(1.0 - costheta*costheta)

# Sample psi.
psi = 2.0*PI*random.random()
cospsi = math.cos(psi)
if (psi < PI):
sinpsi = math.sqrt(1.0 - cospsi*cospsi)
else:
sinpsi = -math.sqrt(1.0 - cospsi*cospsi)

# New trajectory.
if (1 - abs(uz) <= ONE_MINUS_COSZERO) : # close to perpendicular.
uxx = sintheta * cospsi
uyy = sintheta * sinpsi
uzz = costheta * self.SIGN(uz) # SIGN() is faster than division.
else: # usually use this option
temp = math.sqrt(1.0 - uz * uz)
uxx = sintheta * (ux * uz * cospsi - uy * sinpsi) / temp + ux * costheta
uyy = sintheta * (uy * uz * cospsi + ux * sinpsi) / temp + uy * costheta
uzz = -sintheta * cospsi * temp + uz * costheta

# Update trajectory
ux = uxx
uy = uyy
uz = uzz

# Check Roulette
if (W < THRESHOLD):
if (random.random() <= CHANCE):
W = W / CHANCE
else:
photon_status = DEAD
if photon_status is DEAD:
break
if it > max_iterations:
break

if i_photon >= Nphotons:
break
total_reflection = 0.0
for each in range(NR+1):
total_reflection = total_reflection + ReflBin[each]/Nphotons
return total_reflection

def RFresnel(self, n1, n2, cosT1):
r = 0.0
cosT2 = 0.0
COSZERO = 1.0 - 1.0e-12
COS90D = 1 * pow(10,-6)
if n1 == n2: #matched boundary
r = 0.0
cosT2 = cosT1
elif cosT1 > COSZERO: # normal incident
cosT2 = ca1
r = (n2-n1)/(n2+n1)
r *= r
elif cosT1 < COS90D: # very slant
cosT2 = 0.0
r = 1.0
else: #general
sinT1 = math.sqrt(1 - cosT1*cosT1)
sinT2 = n1 * sinT1/n2

if sinT2 >= 1.0:
r = 1.0
cosT2 = 0.0
else:
cosT2 = math.sqrt(1 - sinT2 * sinT2)
cosAP = cosT1*cosT2 - sinT1*sinT2
cosAM = cosT1*cosT2 + sinT1*sinT2
sinAP = sinT1*cosT2 + cosT1*sinT2
sinAM = sinT1*cosT2 - cosT1*sinT2
r = 0.5 * sinAM * sinAM*(cosAM*cosAM+cosAP*cosAP)/(sinAP*sinAP*cosAM*cosAM)
return r

def SIGN(self, x):
if x >=0:
return 1
else:
return 0

skinAbsorption = SkinAbsorption()
reflectanceValues = skinAbsorption.Generate()

import random
import math

import random
import math

class SkinAbsorption:

O2Hb = {}
Hb = {}

def __init__(self):
print(“init SkinAbsorption”)
HbFile = open(‘hb.csv’, “r”)
O2HbFile = open(‘O2Hb.csv’, “r”)

HbLines = HbFile.readlines()
for line in HbLines:
splitLine = line.split(“,”)
self.Hb[int(splitLine[0])] = float(splitLine[1].rstrip(“\n”))

O2HbLines = O2HbFile.readlines()
for line in O2HbLines:
splitLine = line.split(“,”)
self.O2Hb[int(splitLine[0])] = float(splitLine[1].rstrip(“\n”))

HbFile.close()
O2HbFile.close()

def Generate(self):
CSV_Out = “Wavelength,Cm, Ch, Bm,Eumelanin,Pheomelanin,Baseline,Epidermis,ScatteringEpi,ScatteringDerm,Dermis,Reflectance\n”
groupByBlend = []
for Bm in [0.01,0.5,0.99]:
groupByMelanin = []
for Cm in [0.002,0.0135,0.0425,0.1,0.185,0.32,0.5]:
groupByHemoglobin = []
for Ch in [0.003,0.02,0.07,0.16,0.32]:
values = self.GetReflectanceValues( Cm, Ch, Bm )
CSV_Out = CSV_Out + values
groupByHemoglobin.append(values)
groupByMelanin.append(groupByHemoglobin)
groupByBlend.append(groupByMelanin)
CSV_Out_File = open(“AbsorptionValues.csv”,”w+”)
CSV_Out_File.write(CSV_Out)
CSV_Out_File.close()
return groupByBlend

def GetReflectanceValues(self, Cm, Ch, Bm):
CSV_Out = “”
wavelengths = range(380,790,10)
for nm in wavelengths:

# First layer absorption – Epidermis
SAV_eumelanin_L = 6.6 * pow(10,10) * pow(nm,-3.33)
SAV_pheomelanin_L = 2.9 * pow(10,14) * pow(nm,-4.75)
epidermal_hemoglobin_fraction = Ch * 0.25

# baseline – used in both layers
baselineSkinAbsorption_L = 0.0244 + 8.54 * pow(10,-(nm-220)/123)

# Epidermis Absorption Coefficient:
epidermis = Cm * ((Bm * SAV_eumelanin_L) +((1 – Bm ) * SAV_pheomelanin_L)) + ((1 – epidermal_hemoglobin_fraction) * baselineSkinAbsorption_L)

# Second layer absorption – Dermis
gammaOxy = self.O2Hb[int(nm)]
gammaDeoxy = self.Hb[int(nm)]
A = 0.75 * gammaOxy + (1 – 0.75) * gammaDeoxy
dermis = Ch * (A + (( 1 – Ch) * baselineSkinAbsorption_L))

# Scattering coefficients
scattering_epidermis = pow(14.74 * nm, -0.22) + 2.2 * pow(10,11) * pow(nm, -4)
scattering_dermis = scattering_epidermis * 0.5

reflectance = self.MonteCarlo(epidermis, scattering_epidermis, dermis, scattering_dermis, nm)
print(f”Done on {nm} >> ({Cm}, {Ch}, {Bm}) = {reflectance}”)

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”
return CSV_Out

def MonteCarlo (self, epi_mua, epi_mus, derm_mua, derm_mus, nm):
# These are our Monte Carlo Light Transport Variables that don’t change
Nbins = 1000
Nbinsp1 = 1001
PI = 3.1415926
LIGHTSPEED = 2.997925 * pow(10,10)
ALIVE = 1
DEAD = 0
THRESHOLD = 0.01
CHANCE = 0.1
COS90D = 1 * pow(10,-6)
ONE_MINUS_COSZERO = 1 * pow(10,-12)
COSZERO = 1.0 – 1.0e-12 # cosine of about 1e-6 rad
g = 0.9
nt = 1.33 # Index of refraction
epidermis_thickness = 0.25

x = 0.0
y = 0.0
z = 0.0 # photon position

ux = 0.0
uy = 0.0
uz = 0.0 # photon trajectory as cosines

uxx = 0.0
uyy = 0.0
uzz = 0.0 # temporary values used during SPIN

s = 0.0 # step sizes. s = -log(RND)/mus [cm]
costheta = 0.0 # cos(theta)
sintheta = 0.0 # sin(theta)
cospsi = 0.0 # cos(psi)
sinpsi = 0.0 # sin(psi)
psi = 0.0 # azimuthal angle
i_photon = 0.0 # current photon
W = 0.0 # photon weight
absorb = 0.0 # weighted deposited in a step due to absorption
photon_status = 0.0 # flag = ALIVE=1 or DEAD=0
ReflBin = [None] * Nbinsp1 #bin to store weights of escaped photos for reflectivity
epi_albedo = epi_mus/(epi_mus + epi_mua) # albedo of tissue
derm_albedo = derm_mus/(derm_mus + derm_mua) # albedo of tissue
Nphotons = 1000 # number of photons in simulation
NR = Nbins # number of radial positions
radial_size = 2.5 # maximum radial size
r = 0.0 # radial position
dr = radial_size/NR; # cm, radial bin size
ir = 0 # index to radial position
shellvolume = 0.0 # volume of shell at radial position r
CNT = 0.0 # total count of photon weight summed over all bins
rnd = 0.0 # assigned random value 0-1
u = 0.0
temp = 0.0 # dummy variables

# Inits
random.seed(0)
RandomNum = random.random()
for i in range(NR+1):
ReflBin[i] = 0

while True:
i_photon = i_photon + 1

W = 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
while True:
it = it + 1
rnd = random.random()
while rnd <= 0.0: # make sure it is > 0.0
rnd = random.random()
s = -math.log(rnd)/(mua + mus)
x = x + (s * ux)
y = y + (s * uy)
z = z + (s * uz)

if uz < 0:
# calculate partial step to reach boundary surface
s1 = abs(z/uz)
# move back
x = x – (s * ux)
y = y – (s * uy)
z = z – (s * uz)
# take partial step
x = x + (s1 * ux)
y = y + (s1 * uy)
z = z + (s1 * uz)

# photon is now at the surface boundary, figure out how much escaped and how much was reflected
internal_reflectance = self.RFresnel(1.0,nt, -uz )

#Add weighted reflectance of escpaed photon to reflectance bin
external_reflectance = 1 – internal_reflectance
r = math.sqrt(x*x + y*y)
ir = (r/dr)
if ir >= NR:
ir = NR
ReflBin[int(ir)] = ReflBin[int(ir)] + (W * external_reflectance)

# Bounce the photon back into the skin
W = internal_reflectance * W
uz = -uz
x = (s-s1) * ux
y = (s-s1) * uy
z = (s-s1) * uz

# check if we have passed into the second layer, or the first
if z <= epidermis_thickness:
mua = epi_mua
mus = epi_mus
albedo = epi_albedo
else:
mua = derm_mua
mus = derm_mus
albedo = derm_albedo

”’ DROP ”’
absorb = W*(1 – albedo)
W = W – absorb

”’ SPIN ”’
# Sample for costheta
rnd = random.random()
if (g == 0.0):
costheta = 2.0*rnd – 1.0
else:
temp = (1.0 – g*g)/(1.0 – g + 2*g*rnd)
costheta = (1.0 + g*g – temp*temp)/(2.0*g)
sintheta = math.sqrt(1.0 – costheta*costheta)

# Sample psi.
psi = 2.0*PI*random.random()
cospsi = math.cos(psi)
if (psi < PI):
sinpsi = math.sqrt(1.0 – cospsi*cospsi)
else:
sinpsi = -math.sqrt(1.0 – cospsi*cospsi)

# New trajectory.
if (1 – abs(uz) <= ONE_MINUS_COSZERO) : # close to perpendicular.
uxx = sintheta * cospsi
uyy = sintheta * sinpsi
uzz = costheta * self.SIGN(uz) # SIGN() is faster than division.
else: # usually use this option
temp = math.sqrt(1.0 – uz * uz)
uxx = sintheta * (ux * uz * cospsi – uy * sinpsi) / temp + ux * costheta
uyy = sintheta * (uy * uz * cospsi + ux * sinpsi) / temp + uy * costheta
uzz = -sintheta * cospsi * temp + uz * costheta

# Update trajectory
ux = uxx
uy = uyy
uz = uzz

# Check Roulette
if (W < THRESHOLD):
if (random.random() <= CHANCE):
W = W / CHANCE
else:
photon_status = DEAD
if photon_status is DEAD:
break
if it > max_iterations:
break

if i_photon >= Nphotons:
break
total_reflection = 0.0
for each in range(NR+1):
total_reflection = total_reflection + ReflBin[each]/Nphotons
return total_reflection

def RFresnel(self, n1, n2, cosT1):
r = 0.0
cosT2 = 0.0
COSZERO = 1.0 – 1.0e-12
COS90D = 1 * pow(10,-6)
if n1 == n2: #matched boundary
r = 0.0
cosT2 = cosT1
elif cosT1 > COSZERO: # normal incident
cosT2 = ca1
r = (n2-n1)/(n2+n1)
r *= r
elif cosT1 < COS90D: # very slant
cosT2 = 0.0
r = 1.0
else: #general
sinT1 = math.sqrt(1 – cosT1*cosT1)
sinT2 = n1 * sinT1/n2

click to read more

if sinT2 >= 1.0:
r = 1.0
cosT2 = 0.0
else:
cosT2 = math.sqrt(1 – sinT2 * sinT2)
cosAP = cosT1*cosT2 – sinT1*sinT2
cosAM = cosT1*cosT2 + sinT1*sinT2
sinAP = sinT1*cosT2 + cosT1*sinT2
sinAM = sinT1*cosT2 – cosT1*sinT2
r = 0.5 * sinAM * sinAM*(cosAM*cosAM+cosAP*cosAP)/(sinAP*sinAP*cosAM*cosAM)
return r

def SIGN(self, x):
if x >=0:
return 1
else:
return 0

skinAbsorption = SkinAbsorption()
reflectanceValues = skinAbsorption.Generate()

Once its done, create a few charts that span the reflectrance values for 380nm-790nm at varying melanin blend and fractional values. You should see lines that follow these sorts of paths and values:

And that completes our Monte Carlo simulation work. We now need to get this into a usable color.

And that completes our Monte Carlo simulation work. We now need to get this into a usable color.