## Blog Content

Home – Blog Content

# Part 3 (Practical): Spectral Reflection & Color Spaces

## XYZ Spectral Sensitivity Curves

First, let’s get those Gaussian curves for the XYZ spectral sensitivity curves. I’ve pre-generated these into (you guessed it, CSV).

Put this into your development directory, with everything else:

In our __init__ function, do the following:

cmf = open("cie-cmf.csv", "r")
CIE_CMF = [None] * 42
for line in contents:
data = line .rstrip("\n").split(",")
CIE_CMF[int((int(data [0]) - 380) / 10)] = [float(data [1]),float(data [2]),float(data [3])]
cmf.close()
for nm in range(380,790,10):
self.CIE_XYZ_Spectral_Sensitivity_Curve[nm] = CIE_CMF[int((nm - 380)/10)]

And declare at the top of our class, CIE_CMF:

class SkinAbsorption:
O2Hb = {}
Hb = {}
CIE_XYZ_Spectral_Sensitivity_Curve = {}

This will store the values in a dictionary.

## RGB Matrices

Now to add in our RGB matrix and gamma correction for sRGB:

def gamma_correction(self, C):
abs_C = abs(C)
if abs_C > 0.0031308:
return 1.055 * pow(abs_C,1/2.4) - 0.055
else:
return 12.92 * C

def XYZ_to_sRGB(self, xyz):
x = xyz[0]/10
y = xyz[1]/10
z = xyz[2]/10
mat3x3 = [(3.2406, -1.5372, -0.4986), (-0.9689,   1.8758,  0.0415), (0.0557, -0.204,  1.057)]

r = self.gamma_correction(x * mat3x3[0][0] + y * mat3x3[0][1] + z * mat3x3[0][2])
g = self.gamma_correction(x * mat3x3[1][0] + y * mat3x3[1][1] + z * mat3x3[1][2])
b = self.gamma_correction(x * mat3x3[2][0] + y * mat3x3[2][1] + z * mat3x3[2][2])
sRGB = (int(r*255),int(g*255),int(b*255)) #needs to be 0 - 255 for outputing to color image
return sRGB

## Converting Our Data

Our data isn’t in a usable state yet. We need to restructure the data coming from the Monte Carlo simulation to be organised by melanin blend, melanin fraction, hemoglobin fraction – this is so we can output a pixel of colour for each value.

We’ll replace all our CSV output code, and instead of returning a string, we’ll return tuples of the values, like so:

def GetReflectanceValues(self, Cm, Ch, Bm):
wavelengths = range(380,790,10)
reflectances = []
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)
reflectances.append( (reflectance,Cm,Ch,Bm,nm,epidermis,dermis,baselineSkinAbsorption_L))
return reflectances

def Generate(self):
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 )
groupByHemoglobin.append(values)
groupByMelanin.append(groupByHemoglobin)
groupByBlend.append(groupByMelanin)
return groupByBlend

Now, above the return statement add the following to convert all our data to sRGB:

XYZ_Colors = []
specular_col = []
for melanin_blend in groupByBlend:
for melanin_fraction in melanin_blend :
for hemoglobin_fraction in melanin_fraction:
total = (0,0,0)
spec = [0.0] * 41
for data in hemoglobin_fraction:
reflectance = data[0]
spec[int((data[4] -380) / 10)] = data[0]
xyz = self.CIE_XYZ_Spectral_Sensitivity_Curve.get(data[4])
x = xyz[0] * reflectance
y = xyz[1] * reflectance
z = xyz[2] * reflectance
total = (total[0] + x, total[1] + y, total[2] + z)
XYZ_Colors.append(total)

pixelsRGB = []
for xyz in XYZ_Colors:
pixelsRGB.append(self.XYZ_to_sRGB(xyz))

Next step, we’ll make an image!

For reference, current code:

import random
import math
from PIL import Image

class SkinAbsorption:

O2Hb = {}
Hb = {}
CIE_XYZ_Spectral_Sensitivity_Curve = {}

def __init__(self):
print("init SkinAbsorption")

cmf = open("cie-cmf.csv", "r")
CIE_CMF = [None] * 42
for line in contents:
data = line .rstrip("\n").split(",")
CIE_CMF[int((int(data [0]) - 380) / 10)] = [float(data [1]),float(data [2]),float(data [3])]
cmf.close()
for nm in range(380,790,10):
self.CIE_XYZ_Spectral_Sensitivity_Curve[nm] = CIE_CMF[int((nm - 380)/10)]

HbFile = open('hb.csv', "r")
O2HbFile = open('O2Hb.csv', "r")

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

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

HbFile.close()
O2HbFile.close()

def Generate(self):
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 )
groupByHemoglobin.append(values)
groupByMelanin.append(groupByHemoglobin)
groupByBlend.append(groupByMelanin)

XYZ_Colors = []
specular_col = []
for melanin_blend in groupByBlend: # Mb
for melanin_fraction in melanin_blend:
for hemoglobin_fraction in melanin_fraction: # Mf
total = (0,0,0)
spec = [0.0] * 41
for data in hemoglobin_fraction: # Hf

print(data)

reflectance = data[0]

spec[int((data[4] -380) / 10)] = data[0]

xyz = self.CIE_XYZ_Spectral_Sensitivity_Curve.get(data[4])
x = xyz[0] * reflectance
y = xyz[1] * reflectance
z = xyz[2] * reflectance
total = (total[0] + x, total[1] + y, total[2] + z)

XYZ_Colors.append(total)

pixelsRGB = []
for xyz in XYZ_Colors:
pixelsRGB.append(self.XYZ_to_sRGB(xyz))

return pixelsRGB

def GetReflectanceValues(self, Cm, Ch, Bm):
wavelengths = range(380,790,10)
reflectances = []
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)
reflectances.append((reflectance,Cm,Ch,Bm,nm,epidermis,dermis,baselineSkinAbsorption_L))
return reflectances

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
def gamma_correction(self, C):
abs_C = abs(C)
if abs_C > 0.0031308:
return 1.055 * pow(abs_C,1/2.4) - 0.055
else:
return 12.92 * C

def XYZ_to_sRGB(self, xyz):
x = xyz[0]/10
y = xyz[1]/10
z = xyz[2]/10
mat3x3 = [(3.2406, -1.5372, -0.4986), (-0.9689,   1.8758,  0.0415), (0.0557, -0.204,  1.057)]

r = self.gamma_correction(x * mat3x3[0][0] + y * mat3x3[0][1] + z * mat3x3[0][2])
g = self.gamma_correction(x * mat3x3[1][0] + y * mat3x3[1][1] + z * mat3x3[1][2])
b = self.gamma_correction(x * mat3x3[2][0] + y * mat3x3[2][1] + z * mat3x3[2][2])
sRGB = (int(r*255),int(g*255),int(b*255)) #needs to be 0 - 255 for outputing to color image
return sRGB

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

## Most Recent Posts

• All Post
• Articles
• Creative
• Digital
• Marketing
• Quick Tips
• Tangents
• Tutorials
• Uncategorized

2023-04-09

2023-04-09

2023-04-09

2023-04-09