
import os
import shutil
import time
import numpy as np
import matplotlib
from matplotlib import rcParams, style
import matplotlib.pyplot as plt
import pandas as pd


print('numpy version:', np.__version__)
print('matplotlib version:', matplotlib.__version__)
print('pandas version:', pd.__version__)

numpy version: 1.26.4
matplotlib version: 3.9.0
pandas version: 2.2.2


rcParams.update({'font.size': 12})
rcParams.update({'axes.labelsize': 12})
rcParams.update({'axes.titlesize': 12})
rcParams.update({'xtick.labelsize': 12})
rcParams.update({'ytick.labelsize': 12})
rcParams.update({'legend.fontsize': 10})
rcParams.update({'axes.grid': True}) 
plt.rcParams["figure.figsize"] = (6, 4)


import skinoptics
from skinoptics.absorption_coefficient import *
from skinoptics.scattering_coefficient import *
from skinoptics.anisotropy_factor import *
from skinoptics.refractive_index import *
from skinoptics.colors import *


print('pandas version:', skinoptics.__version__)

pandas version: 0.0.1


def mua_EP_Atencio(lambda0, fv_mel):
    return 10*fv_mel*mua_mel_Jacques(lambda0)

def mua_blo_Atencio(lambda0, S):
    return (S*mua_oxy_Prahl(lambda0, Cmass_oxy = 150, molar_mass_oxy = 64500)
            + (1-S)*mua_deo_Prahl(lambda0, Cmass_deo = 150, molar_mass_deo = 64500))

def mua_DE_Atencio(lambda0, fv_blo, S, fv_bil, C_bil):
    return 10*(fv_blo*mua_blo_Atencio(lambda0, S)
               + fv_bil*mua_bil_Li(lambda0, Cmass_bil = C_bil, molar_mass_bil = 585))


all_lambda = np.arange(400, 700+10, 10)
ones = np.ones(len(all_lambda))
color_cycle = ['k', 'r', 'g', 'b', 'c', 'm']


for i in range(3):
    plt.plot(all_lambda, mua_EP_Atencio(all_lambda, 0.02*i),
             color_cycle[i], marker = '.', lw = 2., label = 'fv_mel = {}%'.format(int(0.02*i*100)))
plt.xlabel('wavelength [nm]')
plt.ylabel('absorption coefficient [cm$^{-1}$]')
plt.title('epidermis')
plt.legend(loc = 'upper right')
plt.xlim(400, 700)
plt.ylim(0, 60)
plt.show()


for i in range(6):
    plt.plot(all_lambda, mua_DE_Atencio(all_lambda, 0.002, 0.7, 0.2, 0.05*i),
             color_cycle[i], marker = '.', label = 'C_bil = {:.2f} g/L'.format(int(0.05*i*100)/100))
for i in (1, 3, 5):
    plt.plot(all_lambda, 10*0.2*mua_bil_Li(all_lambda, 0.05*i),
             color_cycle[i], ls = ':', lw = 2., label = 'contribution due to bilirubin only')
plt.xlabel('wavelength [nm]')
plt.ylabel('absorption coefficient [cm$^{-1}$]')
plt.title('dermis (fv_blo = 0.2%, S = 70%, fv_bil = 20%)')
plt.legend(loc = 'upper right')
plt.xlim(400, 700)
plt.ylim(0, 14)
plt.show()


c = 0
for i in [0.001, 0.002, 0.01, 0.02]:
    plt.plot(all_lambda, mua_DE_Atencio(all_lambda, i, 0.7, 0.2, 0.1),
             color_cycle[c], marker = '.', lw = 2., label = 'fv_blo = {} %'.format(i*100))
    c += 1
plt.xlabel('wavelength [nm]')
plt.ylabel('absorption coefficient [cm$^{-1}$]')
plt.title('dermis (C_bil = 0.1 g/L, S = 70%, fv_bil = 20%)')
plt.legend(loc = 'upper right')
plt.xlim(400, 700)
plt.ylim(0, 60)
plt.show()


plt.plot(all_lambda, 10*rmus_Saidi(all_lambda, 6.8, 9.5E10), 
         color_cycle[0], label = 'epidermis & dermis')
plt.plot(all_lambda, 600*ones*(1-0.8), color_cycle[1], lw = 2., label = 'hypodermis')
plt.xlabel('wavelength [nm]')
plt.ylabel('reduced scattering coefficient [cm$^{-1}$]')
plt.legend(loc = 'upper right')
plt.xlim(400, 700)
plt.ylim(0, 160)
plt.show()


plt.plot(all_lambda, 10*rmus_Saidi(all_lambda, 6.8, 9.5E10)/(1 - g_vanGemert(all_lambda)),
         color_cycle[0], label = 'epidermis & dermis')
plt.plot(all_lambda, 600*ones, color_cycle[1], lw = 2., label = 'hypodermis')
plt.xlabel('wavelength [nm]')
plt.ylabel('scattering coefficient [cm$^{-1}$]')
plt.legend(loc = 'upper right')
plt.xlim(400, 700)
plt.ylim(100, 800)
plt.show()


plt.plot(all_lambda, 1.37*ones, color_cycle[0], lw = 2., label = 'epidermis & dermis')
plt.plot(all_lambda, 1.44*ones, color_cycle[1], lw = 2., label = 'hypodermis')
plt.plot(all_lambda, 1.5*ones, color_cycle[2], lw = 2., label = 'skull')
plt.xlabel('wavelength [nm]')
plt.ylabel('refractive index [-]')
plt.legend(loc = 'upper right')
plt.xlim(400, 700)
plt.ylim(1.3, 1.6)
plt.show()


plt.plot(all_lambda, g_vanGemert(all_lambda), color_cycle[0], label = 'epidermis & dermis')
plt.plot(all_lambda, 0.8*ones, color_cycle[1], lw = 2., label = 'hypodermis')
plt.xlabel('wavelength [nm]')
plt.ylabel('anisotropy factor [-]')
plt.legend(loc = 'upper right')
plt.xlim(400, 700)
plt.show()


all_Cbi = [0, 0.05, 0.1, 0.15, 0.20, 0.25]
simulator = 'CUDAMCML' # the user can choose one of the following... 'MCML' or 'CUDAMCML'
folder = r'C:\Users\Lili\Desktop\SkinOptics\simulators\{}'.format(simulator)

# folder2 address
folder2 = r'D:\Victor\2011DelgadoAtencio_{}'.format(simulator)
# create folder2 if it does not exist
if os.path.exists(folder2) == False:
    os.mkdir(folder2)
    
# folder4 address
folder4 = '{}\{}'.format(folder2, 'R_T_A')
# create folder4 if it does not exist
if os.path.exists(folder4) == False:
    os.mkdir(folder4)


'''
This is a python version of the matlab script MCinput.m [DJV11]
'''

for l in range(6):
    Cbi = all_Cbi[l] # Bilirubin content in [g/L]
    
    # 1 SIMULATION PARAMETERS-I ---> The grid and number of photons
    Nf = 50000 # No. of photons
    dz = 20E-4 # dz, dr of the grid
    dr = 20E-4
    Nz = 62 # No. of dz, dr & da
    Nr = 200
    Na = 1

    # 1 SIMULATION PARAMETERS-II---> Geometrical and optical properties of layers
    nc = 3 # Number of layers
    nabove = 1.00 # Refractive index TOP medium
    n1 = 1.37
    d1 = 40E-4 # Layer 1 of 3
    n2 = 1.37
    d2 = 900E-4 # Layer 2 of 3
    n3 = 1.44
    d3 = 300E-4 # Layer 3 of 3
    nbelow = 1.5 # Refractive index BOTTOM medium
    fme = 2/100 # MELANIN content kept constant
    fbl = 0.2/100 # BLOOD content kept constant
    fbi = 20/100 # Partition number
    S = 70/100 # Oxygen saturation

    # 1 SIMULATION PARAMETERS-III---> Spectral data
    lambdai = int(400) # Spectral region and increment
    lambdaf = int(700)
    step = int(10)
    numruns = int(((lambdaf - lambdai)/step) + 1) # Number of runs
    all_lambda = np.arange(lambdai, lambdaf + step, step)  # Discret wavelengths

    # 2 OPTICAL PROPERTIES
    # ABSORPTION and SCATTERING in <<EPIDERMIS>>
    mua1 = mua_EP_Atencio(all_lambda, fv_mel = fme) # Epidermis absorption coefficient
    g = g_vanGemert(all_lambda) # Anisotropy factor
    miuspN = rmus_Saidi(all_lambda, 68, 9.5E11)
    mus1 = miuspN/(1-g) # Epidermis scattering coefficient
    # ABSORPTION and SCATTERING in <<DERMIS>>
    mua2 = mua_DE_Atencio(all_lambda, fv_blo = fbl, S = S, fv_bil = fbi, C_bil = Cbi)
    mus2 = mus1
    # ABSORPTION and SCATTERING in <<HYPODERMIS>>
    mua3 = 1
    mus3 = 600
    g3 = 0.8

    # 3 WRITE INPUT DATA FILES FOR EACH WAVELENGTH 
    # .mci file name
    mci_file = 'Cbi_{}.mci'.format(l)
    # .mci file address
    mci_file_address = '{}/{}'.format(folder, mci_file)
    
    # writting a run for each of the discret wavelengths in the .mci file
    fout = open(mci_file_address, 'w+')
    fout.write('####\n')
    fout.write('# Template of input files for Monte Carlo simulation (mcml).\n')
    fout.write('# Anything in a line after "#" is ignored as comments.\n')
    fout.write('# Space lines are also ignored.\n')
    fout.write('# Lengths are in cm, mua and mus are in 1/cm.\n')
    fout.write('####\n\n')
    fout.write('1.0\t\t\t\t\t\t# file version\n')
    fout.write('{}\t\t\t\t\t\t# number of runs\n\n'.format(numruns))
    for k in np.arange(int(numruns)):
        fout.write('### Specify data for run {}\n'.format(k))
        fout.write('{}.mco\t\tA\t\t\t\t# output filename, ASCII/Binary\n'.format(
            'Cbi_{}_{}'.format(l, all_lambda[k])))
        fout.write('{}\t\t\t\t\t\t# No. of photons\n'.format(Nf))
        fout.write('{}\t{}\t\t\t\t\t# dz, dr\n'.format(dz, dr))
        fout.write('{}\t{}\t{}\t\t\t\t# No. of dz, dr & da.\n\n'.format(Nz, Nr, Na))
        fout.write('{}\t\t\t\t\t\t# No. of layers\n'.format(nc))
        fout.write('# n\tmua\tmus\tg\td\t\t# One line for each layer\n')
        fout.write('{}\t\t\t\t\t\t# n for medium above.\n'.format(nabove))
        fout.write('{}\t{}\t{}\t{}\t{}\t\t# layer 0\n'.format(n1, mua1[k], mus1[k], g[k], d1))
        fout.write('{}\t{}\t{}\t{}\t{}\t\t# layer 1\n'.format(n2, mua2[k], mus2[k], g[k], d2))
        fout.write('{}\t{}\t{}\t{}\t{}\t\t# layer 2\n'.format(n3, mua3, mus3, g3, d3))
        fout.write('{}\t\t\t\t\t\t# n for medium below.\n\n\n'.format(nbelow))
    fout.close()


t0 = time.time()

for l in range(6):
    # part in commom between folder3, .mci and .mco files names
    commom = 'Cbi_{}'.format(l)
    # .mci file name
    mci_file = commom + '.mci'.format(l)
    # .mci file address
    mci_file_address = '{}/{}'.format(folder, mci_file)
    
    # folder3 address
    folder3 = '{}\{}'.format(folder2, commom)
    # create folder3 if it does not exist
    if os.path.exists(folder3) == False:
        os.mkdir(folder3) 
    
    # changing directory to run the Monte Carlo simulator
    os.chdir(folder)
    # running simulator for the .mci file
    os.system('{} {}'.format(simulator, mci_file))
    # returning to previous directory
    os.chdir(r'C:\Users\Lili\Desktop\SkinOptics')
    
    # repeat process for each discret wavelength
    for k in range (lambdai, lambdaf + step, step):
        # .mco file name
        mco_file = '{}_{}.mco'.format(commom, k)
        # .mco file address
        mco_file_address = '{}\{}'.format(folder, mco_file)
        # move .mco file from folder1 to folder3
        shutil.move(mco_file_address, '{}\{}'.format(folder3, mco_file))
    # move .mci file from folder1 to folder2
    shutil.move(mci_file_address, '{}\{}'.format(folder2, mci_file))

print(time.time() - t0)

6.143056154251099


for l in range(6):
    
    # part in commom between folder3 and .mco files names
    commom = 'Cbi_{}'.format(l)
        
    # folder3 address
    folder3 = '{}\{}'.format(folder2, commom)
    # create folder3 if it does not exist
    if os.path.exists(folder3) == False:
        os.mkdir(folder3) 
    
    # .txt file name
    R_T_A_txt_file = 'R_T_A_{}.txt'.format(l)
    # .txt file address
    R_T_A_txt_file_address = '{}\{}'.format(folder4, R_T_A_txt_file)
    # open .txt file
    fout = open(R_T_A_txt_file_address, 'w+')
    # .txt file header
    fout.write('wavelength[nm] Rd[%] T[%] A0[%] A1[%] A2[%] A[%] Rs[%] Rs+Rd+A+T[%] \n')
    
    # repeat process for each discret wavelength
    for k in range (lambdai, lambdaf+step, step):
        
        # .mco file name
        mco_file = '{}_{}.mco'.format(commom, k)
        # .mco file address
        mco_file_address = '{}\{}'.format(folder3, mco_file)
        # open .mco file
        fin = open(mco_file_address, 'r')
        
        # get specular reflectance (Rs) from the .mco file
        for b in range(26):
            line = fin.readline()
        if line.split(' ')[0] == 'RAT': 
            line = fin.readline()
        Rs = float(line.split(' ')[0])*100
        # get diffuse reflectance (Rd) from the .mco file
        line = fin.readline()
        Rd = float(line.split(' ')[0])*100
        # get transmittance (T) from the .mco file
        for b in range(2):
            line = fin.readline()
        T = float(line.split(' ')[0])*100
        for b in range(2):
            line = fin.readline()
        # get absorption (A) from the .mco file
        A = np.zeros(3)
        for b in range(3):
            line = fin.readline()
            A[b] = float(line)*100
        
        # write wavelength, diffuse reflectance, transmittance and absorption
        # and absorption (for each layer and total) to the .txt file
        fout.write('{} {} {} {} {} {} {} {} {}\n'.format(k, Rd, T, A[0], A[1], A[2], A.sum(), Rs, Rd+T+A.sum()+Rs))
        
        # close .mco file
        fin.close()
        
    # close .txt file
    fout.close()


Rd_MCML, Rd_CUDAMCML = np.zeros([2, 6, len(all_lambda)])

for i in range(6):
    Rd_MCML[i,:] = np.array(pd.read_csv(r'D:\Victor\2011DelgadoAtencio_{}\R_T_A\R_T_A_{}.txt'.format('MCML', i), sep = ' '))[:,1]
    Rd_CUDAMCML[i,:] = np.array(pd.read_csv(r'D:\Victor\2011DelgadoAtencio_{}\R_T_A\R_T_A_{}.txt'.format('CUDAMCML', i), sep = ' '))[:,1]


i = 0
plt.scatter(all_lambda, Rd_MCML[i,:], color = color_cycle[i], marker = '.', label = 'MCML')
for i in range(6):
    plt.plot(all_lambda, Rd_CUDAMCML[i,:], color = color_cycle[i], lw = 2., label = 'CUDAMCML, C_bil = {:.2f} g/L'.format(all_Cbi[i]))
    plt.scatter(all_lambda, Rd_MCML[i,:], color = color_cycle[i], marker = '.')
plt.xlabel('wavelength [nm]')
plt.ylabel('diffuse reflectance [%]')
plt.title('(fv_mel = 2%, fv_blo = 0.2%, S = 70%, fv_bil = 20%)')
plt.legend(loc = 'lower right')
plt.xlim(400, 700)
plt.ylim(0, 70)
plt.show()


for i in range(6):
    plt.plot(all_lambda, np.abs(Rd_CUDAMCML[i,:] - Rd_MCML[i,:]),
             color = color_cycle[i], marker = '.', lw = 2., label = 'C_bil = {:.2f} g/L'.format(all_Cbi[i]))
plt.xlabel('wavelength [nm]')
plt.ylabel('absolute difference in Rd [%]')
plt.title('CUDAMCML & MCML')
plt.legend()
plt.xlim(400, 700)
plt.ylim(0, 2.5)
plt.show()


T_MCML, T_CUDAMCML = np.zeros([2, 6, len(all_lambda)])

for i in range(6):
    T_MCML[i,:] = np.array(pd.read_csv(r'D:\Victor\2011DelgadoAtencio_{}\R_T_A\R_T_A_{}.txt'.format('MCML', i), sep = ' '))[:,2]
    T_CUDAMCML[i,:] = np.array(pd.read_csv(r'D:\Victor\2011DelgadoAtencio_{}\R_T_A\R_T_A_{}.txt'.format('CUDAMCML', i), sep = ' '))[:,2]


i = 0
plt.scatter(all_lambda, T_MCML[i,:], color = color_cycle[i], marker = '.', label = 'MCML')
for i in range(6):
    plt.plot(all_lambda, T_CUDAMCML[i,:],
             color = color_cycle[i], lw = 2., label = 'CUDAMCML, C_bil = {:.2f} g/L'.format(all_Cbi[i]))
    plt.scatter(all_lambda, T_MCML[i,:], color = color_cycle[i], marker = '.')
plt.xlabel('wavelength [nm]')
plt.ylabel('transmittance [%]')
plt.title('(fv_mel = 2%, f_blo = 0.2%, S = 70%, fv_bil = 20%)')
plt.legend(loc = 'upper left')
plt.xlim(400, 700)
plt.ylim(0, 50)
plt.show()


for i in range(6):
    plt.plot(all_lambda, np.abs(T_CUDAMCML[i,:] - T_MCML[i,:]),
                color = color_cycle[i], marker = '.', lw = 2., label = 'C_bil = {:.2f} g/L'.format(all_Cbi[i]))
plt.xlabel('wavelength [nm]')
plt.ylabel('absolute difference in T [%]')
plt.title('CUDAMCML & MCML')
plt.legend()
plt.xlim(400, 700)
plt.ylim(0, 2.5)
plt.show()


RTA_sum_MCML, RTA_sum_CUDAMCML = np.zeros([2, 6, len(all_lambda)])

for i in range(6):
    RTA_sum_MCML[i,:] = np.array(pd.read_csv(r'D:\Victor\2011DelgadoAtencio_{}\R_T_A\R_T_A_{}.txt'.format('MCML', i), sep = ' '))[:,8]
    RTA_sum_CUDAMCML[i,:] = np.array(pd.read_csv(r'D:\Victor\2011DelgadoAtencio_{}\R_T_A\R_T_A_{}.txt'.format('CUDAMCML', i), sep = ' '))[:,8]


i = 0
plt.scatter(all_lambda, RTA_sum_MCML[i,:],
            color = color_cycle[i], marker = '.', label = 'MCML')
for i in range(6):
    plt.plot(all_lambda, RTA_sum_CUDAMCML[i,:],
             color = color_cycle[i], lw = 2., label = 'CUDAMCML, C_bil = {:.2f} g/L'.format(all_Cbi[i]))
    plt.scatter(all_lambda, RTA_sum_MCML[i,:],
                color = color_cycle[i], marker = '.')
plt.ylabel('R+T+A [%]')
plt.title('(fv_mel = 2%, f_blo = 0.2%, S = 70%, fv_bil = 20%)')
plt.xlim(400, 700)
plt.ylim(98.9, 100.1)
plt.legend(loc = 'lower right')
plt.show()


pd.DataFrame(np.array([[RTA_sum_MCML.min(), RTA_sum_MCML.max()],
                       [RTA_sum_CUDAMCML.min(), RTA_sum_CUDAMCML.max()]]),
             index = ['R+T+A (MCML) [%]', 'R+T+A (CUDAMCML) [%]'],
             columns = ['min', 'max'])


R, G, B = np.zeros([3, 6])

for i in range(6):
    R[i], G[i], B[i] = sRGB_from_spectrum(all_lambda, Rd_CUDAMCML[i,:])
RGB = np.vstack((R, G, B)).T


rcParams.update({'axes.grid': False})
for i in range(6):
    plt.gca().add_patch(plt.Rectangle((i, 0.), 1., 1., 
                                      fc = RGB[i], 
                                      ec = (0, 0, 0)))
plt.xticks([0.5, 1.5, 2.5, 3.5, 4.5, 5.5], all_Cbi)
plt.yticks([0], '')
plt.xlabel('bilirubin concentration [g/L]')
plt.axis('scaled')
plt.axis([0, 6, 0, 1])
plt.show()
rcParams.update({'axes.grid': True})

Reproducing results obtained by Delgado Atencio, Jacques & Vázquez y Montiel 2011 in Python with MCML and CUDAMCML¶

Absorption Coefficient¶

Scattering Coefficient¶

Refractive Index¶

Anisotropy Factor¶

Simulations¶

Diffuse Reflectance¶

Extra: Transmittance¶

Extra: Checking if R+T+A = 100%¶

Extra: sRGB color space¶

	min	max
R+T+A (MCML) [%]	99.99129	100.008290
R+T+A (CUDAMCML) [%]	99.50065	99.924765