## half4.xyz  ## Blog Content

Home – Blog Content

# 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: σepidermisa, σepidermiss, σdermisa, σdermiss 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``````
``````def SIGN(self, x):
if x >=0:
return 1
else:
return 0  ``````

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
THRESHOLD = 0.01
CHANCE = 0.1
COS90D = exp(1,-6)
ONE_MINUS_COSZERO = exp(1,-12)
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
r = 0.0 # radial position
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

``````

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:
break
if it > max_iterations:
break``````

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

``````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)``````

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``````

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``````

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

``````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``````

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

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:

``````[....]
# 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()

``````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.

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")

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

for line in O2HbLines:
splitLine = line.split(",")
self.O2Hb[int(splitLine)] = float(splitLine.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)]
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
THRESHOLD = 0.01
CHANCE = 0.1
COS90D = 1 * pow(10,-6)
ONE_MINUS_COSZERO = 1 * pow(10,-12)
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
r = 0.0 # radial position
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:
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

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()``````

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.

## Most Recent Posts

• All Post
• Articles
• Creative
• Digital
• Marketing
• Quick Tips
• Tangents
• Tutorials
• Uncategorized 2023-04-09 2023-04-09 ## Digital Marketing Made Easy: Let Our Team Handle

2023-04-09 2023-04-09