![\"\"](\"https:\/\/half4.xyz\/wp-content\/uploads\/2020\/07\/image-38.png\")
<\/figure>\n\n\n\nIn our __init__ function, do the following:<\/p>\n\n\n\n
cmf = open(\"cie-cmf.csv\", \"r\") \ncontents = cmf.readlines()\nCIE_CMF = [None] * 42\nfor line in contents:\n data = line .rstrip(\"\\n\").split(\",\")\n CIE_CMF[int((int(data [0]) - 380) \/ 10)] = [float(data [1]),float(data [2]),float(data [3])]\ncmf.close()\nfor nm in range(380,790,10):\n self.CIE_XYZ_Spectral_Sensitivity_Curve[nm] = CIE_CMF[int((nm - 380)\/10)]<\/code><\/pre>\n\n\n\n<\/p>\n\n\n\n
And declare at the top of our class, CIE_CMF:<\/p>\n\n\n\n
class SkinAbsorption:\n O2Hb = {}\n Hb = {}\n CIE_XYZ_Spectral_Sensitivity_Curve = {}<\/code><\/pre>\n\n\n\n
CIE_XYZ_Spectral_Sensitivity_Curve = {}<\/code><\/p>\n\n\n\nThis will store the values in a dictionary.<\/p>\n\n\n\n
RGB Matrices<\/h2>\n\n\n\n
Now to add in our RGB matrix and gamma correction for sRGB:<\/p>\n\n\n\n
def gamma_correction(self, C):\n abs_C = abs(C)\n if abs_C > 0.0031308:\n return 1.055 * pow(abs_C,1\/2.4) - 0.055\n else:\n return 12.92 * C\n\ndef XYZ_to_sRGB(self, xyz):\n x = xyz[0]\/10\n y = xyz[1]\/10\n z = xyz[2]\/10\n mat3x3 = [(3.2406, -1.5372, -0.4986), (-0.9689, 1.8758, 0.0415), (0.0557, -0.204, 1.057)]\n \n r = self.gamma_correction(x * mat3x3[0][0] + y * mat3x3[0][1] + z * mat3x3[0][2])\n g = self.gamma_correction(x * mat3x3[1][0] + y * mat3x3[1][1] + z * mat3x3[1][2])\n b = self.gamma_correction(x * mat3x3[2][0] + y * mat3x3[2][1] + z * mat3x3[2][2])\n sRGB = (int(r*255),int(g*255),int(b*255)) #needs to be 0 - 255 for outputing to color image\n return sRGB<\/code><\/pre>\n\n\n\n<\/p>\n\n\n\n
Converting Our Data<\/h2>\n\n\n\n
Our data isn\u2019t 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 \u2013 this is so we can output a pixel of colour for each value.<\/p>\n\n\n\n
We\u2019ll replace all our CSV output code, and instead of returning a string, we\u2019ll return tuples of the values, like so:<\/p>\n\n\n\n
def GetReflectanceValues(self, Cm, Ch, Bm):\n wavelengths = range(380,790,10)\n reflectances = []\n for nm in wavelengths:\n\t\t\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\t\t\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\t\t\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\t\t\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\t\t\n reflectance = self.MonteCarlo(epidermis, scattering_epidermis, dermis, scattering_dermis, nm)\n reflectances.append( (reflectance,Cm,Ch,Bm,nm,epidermis,dermis,baselineSkinAbsorption_L))\nreturn reflectances<\/code><\/pre>\n\n\n\n <\/code><\/p>\n\n\n\nAnd adjust the generate function:<\/p>\n\n\n\n
def Generate(self):\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 groupByHemoglobin.append(values)\n groupByMelanin.append(groupByHemoglobin)\n groupByBlend.append(groupByMelanin)\n return groupByBlend<\/code><\/pre>\n\n\n\n<\/p>\n\n\n\n
Now, above the return statement add the following to convert all our data to sRGB:<\/p>\n\n\n\n
XYZ_Colors = []\nspecular_col = []\nfor melanin_blend in groupByBlend:\n for melanin_fraction in melanin_blend :\n for hemoglobin_fraction in melanin_fraction:\n total = (0,0,0)\n spec = [0.0] * 41\n for data in hemoglobin_fraction:\n reflectance = data[0]\n spec[int((data[4] -380) \/ 10)] = data[0]\n xyz = self.CIE_XYZ_Spectral_Sensitivity_Curve.get(data[4])\n x = xyz[0] * reflectance\n y = xyz[1] * reflectance\n z = xyz[2] * reflectance\n total = (total[0] + x, total[1] + y, total[2] + z)\n XYZ_Colors.append(total)\n\n\npixelsRGB = []\nfor xyz in XYZ_Colors:\n pixelsRGB.append(self.XYZ_to_sRGB(xyz))<\/code><\/pre>\n\n\n\n<\/p>\n\n\n\n
Next step, we\u2019ll make an image! <\/p>\n\n\n\n
For reference, current code:<\/p>\n\n\n\n
import random\nimport math\nfrom PIL import Image\n\nclass SkinAbsorption:\n\n O2Hb = {}\n Hb = {}\n CIE_XYZ_Spectral_Sensitivity_Curve = {}\n \n def __init__(self):\n print(\"init SkinAbsorption\")\n\n cmf = open(\"cie-cmf.csv\", \"r\") \n contents = cmf.readlines()\n CIE_CMF = [None] * 42\n for line in contents:\n data = line .rstrip(\"\\n\").split(\",\")\n CIE_CMF[int((int(data [0]) - 380) \/ 10)] = [float(data [1]),float(data [2]),float(data [3])]\n cmf.close()\n for nm in range(380,790,10):\n self.CIE_XYZ_Spectral_Sensitivity_Curve[nm] = CIE_CMF[int((nm - 380)\/10)]\n \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 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 groupByHemoglobin.append(values)\n groupByMelanin.append(groupByHemoglobin)\n groupByBlend.append(groupByMelanin)\n\n XYZ_Colors = []\n specular_col = []\n for melanin_blend in groupByBlend: # Mb\n for melanin_fraction in melanin_blend:\n for hemoglobin_fraction in melanin_fraction: # Mf\n total = (0,0,0)\n spec = [0.0] * 41\n for data in hemoglobin_fraction: # Hf\n \n print(data)\n \n reflectance = data[0]\n\n spec[int((data[4] -380) \/ 10)] = data[0]\n\n xyz = self.CIE_XYZ_Spectral_Sensitivity_Curve.get(data[4])\n x = xyz[0] * reflectance\n y = xyz[1] * reflectance\n z = xyz[2] * reflectance\n total = (total[0] + x, total[1] + y, total[2] + z)\n\n XYZ_Colors.append(total)\n\n\n pixelsRGB = []\n for xyz in XYZ_Colors:\n pixelsRGB.append(self.XYZ_to_sRGB(xyz))\n\n return pixelsRGB\n \n def GetReflectanceValues(self, Cm, Ch, Bm):\n wavelengths = range(380,790,10)\n reflectances = []\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 reflectances.append((reflectance,Cm,Ch,Bm,nm,epidermis,dermis,baselineSkinAbsorption_L))\n return reflectances\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 def gamma_correction(self, C):\n abs_C = abs(C)\n if abs_C > 0.0031308:\n return 1.055 * pow(abs_C,1\/2.4) - 0.055\n else:\n return 12.92 * C\n\n def XYZ_to_sRGB(self, xyz):\n x = xyz[0]\/10\n y = xyz[1]\/10\n z = xyz[2]\/10\n mat3x3 = [(3.2406, -1.5372, -0.4986), (-0.9689, 1.8758, 0.0415), (0.0557, -0.204, 1.057)]\n \n r = self.gamma_correction(x * mat3x3[0][0] + y * mat3x3[0][1] + z * mat3x3[0][2])\n g = self.gamma_correction(x * mat3x3[1][0] + y * mat3x3[1][1] + z * mat3x3[1][2])\n b = self.gamma_correction(x * mat3x3[2][0] + y * mat3x3[2][1] + z * mat3x3[2][2])\n sRGB = (int(r*255),int(g*255),int(b*255)) #needs to be 0 - 255 for outputing to color image\n return sRGB\n\nskinAbsorption = SkinAbsorption()\nreflectanceValues = skinAbsorption.Generate()<\/code><\/pre>\n\n\n\n<\/p>\n","protected":false},"excerpt":{"rendered":"
XYZ Spectral Sensitivity Curves First, let\u2019s get those Gaussian curves for the XYZ spectral sensitivity curves. I\u2019ve pre-generated these into (you guessed it, CSV). Put this into your development directory, with everything else: In our __init__ function, do the following: And declare at the top of our class, CIE_CMF: 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: Converting Our Data Our data isn\u2019t 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 \u2013 this is so we can output a pixel of colour for each value. We\u2019ll replace all our CSV output code, and instead of returning a string, we\u2019ll return tuples of the values, like so: And adjust the generate function: Now, above the return statement add the following to convert all our data to sRGB: Next step, we\u2019ll make an image! For reference, current code:<\/p>\n","protected":false},"author":1,"featured_media":297,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[1],"tags":[8,7],"class_list":["post-152","post","type-post","status-publish","format-standard","has-post-thumbnail","hentry","category-uncategorized","tag-lut","tag-skin"],"_links":{"self":[{"href":"https:\/\/half4.xyz\/index.php\/wp-json\/wp\/v2\/posts\/152","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/half4.xyz\/index.php\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/half4.xyz\/index.php\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/half4.xyz\/index.php\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/half4.xyz\/index.php\/wp-json\/wp\/v2\/comments?post=152"}],"version-history":[{"count":14,"href":"https:\/\/half4.xyz\/index.php\/wp-json\/wp\/v2\/posts\/152\/revisions"}],"predecessor-version":[{"id":1300,"href":"https:\/\/half4.xyz\/index.php\/wp-json\/wp\/v2\/posts\/152\/revisions\/1300"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/half4.xyz\/index.php\/wp-json\/wp\/v2\/media\/297"}],"wp:attachment":[{"href":"https:\/\/half4.xyz\/index.php\/wp-json\/wp\/v2\/media?parent=152"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/half4.xyz\/index.php\/wp-json\/wp\/v2\/categories?post=152"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/half4.xyz\/index.php\/wp-json\/wp\/v2\/tags?post=152"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}