diff --git a/Modules/Biophotonics/python/generateNoisyRandomSpectra.py b/Modules/Biophotonics/python/generateNoisyRandomSpectra.py index 679b1370d8..8f48790d92 100644 --- a/Modules/Biophotonics/python/generateNoisyRandomSpectra.py +++ b/Modules/Biophotonics/python/generateNoisyRandomSpectra.py @@ -1,88 +1,88 @@ # -*- coding: utf-8 -*- """ Created on Fri Feb 6 16:14:57 2015 @author: wirkert """ import time import datetime import random import numpy as np import monteCarloHelper as mch import setupSimulation as setup - # the input file without the run specific parameters for ua, us and d: -infileString = 'input/darmpythontest.mci' +infileString = 'data/colonTemplate.mci' infile = open(infileString) # the output folder for the mc simulations -outfolderMC = 'outputMC/' +# attention: this is relative to your gpumcml path! +outfolderMC ='outputMC/' # the output folder for the reflectance spectra -outfolderRS = 'outputRS/' - +outfolderRS = 'data/output/' +gpumcmlDirectory = '/home/wirkert/workspace/monteCarlo/gpumcml/fast-gpumcml/' +gpumcmlExecutable = 'gpumcml.sm_20' BVFs, Vss, ds, SaO2s, rs, nrSamples, photons, wavelengths, FWHM, eHbO2, eHb = setup.setupNormalSimulation() nrSimulations = 100 reflectances = np.zeros((nrSimulations, len(wavelengths))) parameters = np.zeros((nrSimulations, 8)) print('start simulations...') #%% start program logic start = time.time() for i in range(nrSimulations): - j = 0 + print('starting simulation ' + str(i) + ' of ' + str(nrSimulations)) BVF = random.uniform(min(BVFs), max(BVFs)) Vs = random.uniform(min(Vss), max(Vss)) d = random.uniform(min(ds), max(ds)) r = random.uniform(min(rs), max(rs)) SaO2= random.uniform(min(SaO2s), max(SaO2s)) sm_BVF = random.uniform(min(BVFs), max(BVFs)) sm_Vs = random.uniform(min(Vss), max(Vss)) sm_SaO2= random.uniform(min(SaO2s), max(SaO2s)) parameters[i,:] = np.array([BVF, Vs, d, r, SaO2, sm_BVF, sm_Vs, sm_SaO2]) - for wavelength in wavelengths: + for j, wavelength in enumerate(wavelengths): reflectanceValue = mch.runOneSimulation( wavelength, eHbO2, eHb, - infile, outfolderMC, + infile, outfolderMC, gpumcmlDirectory, gpumcmlExecutable, BVF, Vs, d, r, SaO2, submucosa_BVF=sm_BVF, submucosa_Vs=sm_Vs, submucosa_SaO2=sm_SaO2, Fwhm = FWHM, nrPhotons=photons) print((BVF, Vs, d, SaO2, r, wavelength, sm_BVF, sm_Vs, sm_SaO2)) # here, summarize result from wavelength in reflectance spectrum reflectances[i, j] = reflectanceValue - j +=1 infile.close() # save the reflectance results! now = datetime.datetime.now().strftime("%Y%B%d%I:%M%p") np.save(outfolderRS + now + "reflectancesRandomWithNoise", reflectances) np.save(outfolderRS + now + "paramterersRandomWithNoise", parameters) end = time.time() print "total time for script: " + str((end - start)) \ No newline at end of file diff --git a/Modules/Biophotonics/python/generateRandomSpectra.py b/Modules/Biophotonics/python/generateRandomSpectra.py index f9dd66db28..98b5382421 100644 --- a/Modules/Biophotonics/python/generateRandomSpectra.py +++ b/Modules/Biophotonics/python/generateRandomSpectra.py @@ -1,82 +1,93 @@ # -*- coding: utf-8 -*- """ Created on Fri Feb 6 16:14:57 2015 @author: wirkert """ import time import datetime import random import numpy as np import monteCarloHelper as mch import setupSimulation as setup # the input file without the run specific parameters for ua, us and d: infileString = 'input/darmpythontest.mci' infile = open(infileString) # the output folder for the mc simulations outfolderMC = 'outputMC/' # the output folder for the reflectance spectra outfolderRS = 'outputRS/' +# the input file without the run specific parameters for ua, us and d: +infileString = 'data/colonTemplate.mci' +infile = open(infileString) +# the output folder for the mc simulations +# attention: this is relative to your gpumcml path! +outfolderMC ='outputMC/' +# the output folder for the reflectance spectra +outfolderRS = 'data/output/' +gpumcmlDirectory = '/home/wirkert/workspace/monteCarlo/gpumcml/fast-gpumcml/' +gpumcmlExecutable = 'gpumcml.sm_20' + + BVFs, Vss, ds, SaO2s, rs, nrSamples, photons, wavelengths, FWHM, eHbO2, eHb = setup.setupNormalSimulation() nrSimulations = 100 reflectances = np.zeros((nrSimulations, len(wavelengths))) parameters = np.zeros((nrSimulations, 2)) print('start simulations...') #%% start program logic start = time.time() for i in range(nrSimulations): - j = 0 + print('starting simulation ' + str(i) + ' of ' + str(nrSimulations)) BVF = random.uniform(min(BVFs), max(BVFs)) Vs = random.uniform(min(Vss), max(Vss)) d = random.uniform(min(ds), max(ds)) r = random.uniform(min(rs), max(rs)) SaO2= random.uniform(min(SaO2s), max(SaO2s)) parameters[i,:] = np.array([BVF, Vs])#, d, r, SaO2]) - for wavelength in wavelengths: + for j, wavelength in enumerate(wavelengths): reflectanceValue = mch.runOneSimulation( wavelength, eHbO2, eHb, - infile, outfolderMC, + infile, outfolderMC, gpumcmlDirectory, gpumcmlExecutable, BVF, Vs, d, r, SaO2, Fwhm = FWHM, nrPhotons=photons) print((BVF, Vs, d, SaO2, r, wavelength)) # here, summarize result from wavelength in reflectance spectrum reflectances[i, j] = reflectanceValue - j +=1 infile.close() # save the reflectance results! now = datetime.datetime.now().strftime("%Y%B%d%I:%M%p") np.save(outfolderRS + now + "reflectancesRandom", reflectances) np.save(outfolderRS + now + "paramtersRandom", parameters) end = time.time() print "total time for script: " + str((end - start)) \ No newline at end of file diff --git a/Modules/Biophotonics/python/monteCarloHelper.py b/Modules/Biophotonics/python/monteCarloHelper.py index 96e728d870..1e0975d74e 100644 --- a/Modules/Biophotonics/python/monteCarloHelper.py +++ b/Modules/Biophotonics/python/monteCarloHelper.py @@ -1,160 +1,160 @@ # -*- coding: utf-8 -*- """ Created on Sun Feb 01 13:52:23 2015 Some helper methods @author: Seb """ import math import numpy as np import mieMonteCarlo import subprocess import os import contextlib @contextlib.contextmanager def cd(newPath): savedPath = os.getcwd() os.chdir(newPath) yield os.chdir(savedPath) def createSimulationFile(infile, replacements, simulationFileName): infile.seek(0); simulationFile = open(simulationFileName, 'w') for line in infile: for src, target in replacements.iteritems(): line = line.replace(src, target) simulationFile.write(line) simulationFile.close() def calcUa(BVF, cHb, SaO2, eHbO2, eHb, wavelength): # determine ua [1/m] as combination of extinction coefficients. # For more on this equation, please refer to # http://omlc.org/spectra/hemoglobin/ return BVF * math.log(10) * cHb * \ (SaO2 * eHbO2(wavelength) + (1-SaO2) * eHb(wavelength)) \ / 64500 def calcUsMuscle(wavelength): return 168.52 * (wavelength * 10**9)**-0.332 / (1 - 0.96) * 100 def getReflectanceFromFile(fileName): ''' extract reflectance from a specific monte carlo output file. Attention: this is kept very simple and assumes that the reflectance output starts at a specific line and keeps till eof. This might be subject to change, e.g., when you alter the numbers of layers the line number there reflectance values start will change. Error handling etc. needed. ''' with open(fileName) as myFile: for num, line in enumerate(myFile, 1): if "Rd_r #Rd[0], [1],..Rd[nr-1]. [1/cm2]" in line: break #print 'found at line:', num reflectances = np.loadtxt(fileName, skiprows=num) return sum(reflectances) def normalizeIntegral(data, wavelengths): # normalize data norms = np.trapz(data, wavelengths, axis=1) return data / norms[:,None] def normalizeImageQuotient(data, iqBand=0): # use line 0 as image quotient normData = data / data[:,iqBand][:,None] # discard it normData = np.concatenate((normData[:,0:iqBand], normData[:,iqBand+1:]), axis=1) return normData def runOneSimulation(wavelength, eHbO2, eHb, - infile, simulationOutputFolderName, gpumcmlDirectory, gpymcmlExecutable, + infile, simulationOutputFolderName, gpumcmlDirectory, gpumcmlExecutable, BVF=0.1, Vs=0.4, d=500 * 10**-6, r=0.4 * 10**-6, SaO2=0.7, cHb=120, submucosa_BVF=0.1, submucosa_Vs=0.3, submucosa_SaO2=0.7, Fwhm = 20 * 10**-9, nrPhotons = 10**6, ): """ Args: _____ gpumcmlString: a string containing the full absolute path gpumcml (e.g. /home/wirkert/gpumcml.sm20) cHb: concentration of haemoglobin per unit volume of blood in colon [g/L], taken from: "Modelling and validation of spectral reflectance for the colon" """ outputFilename = \ simulationOutputFolderName + \ "BVF" + repr(BVF) + \ "Vs" + repr(Vs) + \ "d" + repr(d) + \ "SaO2" + repr(SaO2) + \ "r" + repr(r) + \ "photons" + repr(nrPhotons) + \ "wavelength" + repr(wavelength) + \ ".mco" usg = mieMonteCarlo.mieMonteCarlo_FWHM(wavelength, r, Vs, FWHM = Fwhm) ua = calcUa(BVF, cHb, SaO2, eHbO2, eHb, wavelength) # calculate ua, us for submucosa, values taken from # "Model based inversion for deriving maps of # histological parameters characteristic of cancer # from ex-vivo multispectral images of the colon" ua_sm = calcUa(submucosa_BVF, cHb, submucosa_SaO2, eHbO2, eHb, wavelength) usg_sm = mieMonteCarlo.mieMonteCarlo_FWHM(wavelength, 2*10**-6, submucosa_Vs, FWHM = Fwhm) # now create us and ua for muscle layer according to # "Modelling and validation of spectral reflectance for the colon" us_mus = calcUsMuscle(wavelength) A = 1.7923385088285804 # calculated to retrieve an absorption of 11.2 cm-1 # at 515nm ua_mus = calcUa(A*0.1, 120, 0.7, eHbO2, eHb, wavelength) # factors are due to the program expecting [cm] instead of [m] replacements = {'muc_c_ua':str(ua / 100), 'muc_c_us':str(usg['us'] / 100), 'muc_c_g':str(usg['g']), 'muc_c_d':str(d * 100), 'sm_c_ua':str(ua_sm / 100), 'sm_c_us':str(usg_sm['us'] / 100), 'sm_c_g':str(usg_sm['g']), 'mus_c_ua':str(ua_mus / 100), 'mus_c_us':str(us_mus / 100), 'outputFilename':outputFilename, 'c_ph': str(nrPhotons)} createSimulationFile(infile, replacements, "data/output/simulationFile.mci") - args = ("./" + gpymcmlExecutable, "-A", os.getcwd() + "/data/output/simulationFile.mci") + args = ("./" + gpumcmlExecutable, "-A", os.getcwd() + "/data/output/simulationFile.mci") with cd(gpumcmlDirectory): popen = subprocess.Popen(args, stdout=subprocess.PIPE) popen.wait() # outside the context manager we are back wherever we started. # this function is error prone, please refer to documentation to see # if problematic for you return getReflectanceFromFile(gpumcmlDirectory + outputFilename) if __name__ == "__main__": us515 = calcUsMuscle(515 * 10**-9) print ("calculated us [cm-1] for muscle layer at 515nm: " + str(us515 / 100) + " expected ca. 530") us1064 = calcUsMuscle(1064 * 10**-9) print ("calculated us [cm-1] for muscle layer at 1064nm: " + str(us1064 / 100)) \ No newline at end of file diff --git a/Modules/Biophotonics/python/optimization.py b/Modules/Biophotonics/python/optimization.py index 8834c5d26d..dffd7e2816 100644 --- a/Modules/Biophotonics/python/optimization.py +++ b/Modules/Biophotonics/python/optimization.py @@ -1,162 +1,148 @@ # -*- coding: utf-8 -*- """ Created on Sat Feb 07 23:42:41 2015 @author: Seb """ from scipy.interpolate import RectBivariateSpline import numpy as np import matplotlib.pyplot as plt from scipy.optimize import minimize from reflectanceError import ReflectanceError from mpl_toolkits.mplot3d import Axes3D from matplotlib import cm -import setupSimulation as setup - - - -#BVFs, Vss, ds, SaO2s, rs, nrSamples, photons, wavelengths, FWHM, eHbO2, eHb = setup.setupNormalSimulation() - - - - dataFolder = "outputRS/" # load data trainingParameters = np.load(dataFolder + "paramterers2015February0511:02PM.npy") trainingReflectances = np.load(dataFolder + "reflectances2015February0511:02PM.npy") BVFs = np.unique(trainingParameters[:,0]) Vss = np.unique(trainingParameters[:,1]) reflectanceGrid3D = np.reshape(trainingReflectances, (len(Vss), len(BVFs), trainingReflectances.shape[1]), order='F') functionToMinimize = ReflectanceError(Vss, BVFs, reflectanceGrid3D) testParameters = np.load(dataFolder + "2015February1107:43PMparamterersRandom.npy") testReflectances = np.load(dataFolder + "2015February1107:43PMreflectancesRandom.npy") absErrors = np.zeros_like(testParameters) functionToMinimize.setReflectanceToMatch(trainingReflectances[0,:]) absError = minimize(functionToMinimize.f, [0.3, 0.05], method="Nelder-Mead") #absError = minimize(functionToMinimize.f, [0.04, 0.3], method="BFGS", jac=functionToMinimize.df) print(" function evaluation yields: " + str(absError.x) + " with absolute error: " + str(functionToMinimize.f(absError.x))) print(" expected evaluation result: " + str([0.04, 0.01]) + " real result error: " + str(functionToMinimize.f([0.04, 0.01]))) for idx, (testParameter, testReflectance) in enumerate(zip(testParameters, testReflectances)): functionToMinimize.setReflectanceToMatch(testReflectance) minimization = minimize(functionToMinimize.f, [0.3, 0.05], method="Nelder-Mead") # interpolation extrapolates with constant values. Since the minimization function is # covex, we can just crop it to the bounds # todo: sort out this ugly mess of flipping Vss and BVFs! estimatedVs = np.clip(minimization.x[0], min(Vss), max(Vss)) estimatedBVF= np.clip(minimization.x[1], min(BVFs), max(BVFs)) - print (str(idx) + ": estimated: " + str([[estimatedBVF, estimatedVs]]) + " expected: " + str(testParameter)) - print (str(testReflectance)) - print (" match error: " + str(functionToMinimize.f(minimization.x)) + \ - " real error: " + str(functionToMinimize.f([testParameter[1], testParameter[0]]))) - absErrors[idx,:] = np.abs([estimatedBVF, estimatedVs] - testParameter) # plot print("error distribution BVF, Volume fraction") print("median: " + str(np.median(absErrors, axis=0))) print("lower quartile: " + str(np.percentile(absErrors, 25, axis=0))) print("higher quartile: " + str(np.percentile(absErrors, 75, axis=0))) if __name__ == "__main__": rbs = RectBivariateSpline(Vss, BVFs, np.reshape(trainingReflectances[:,0],(100,100)).T) rbsFrom3DGrid = RectBivariateSpline(Vss, BVFs, reflectanceGrid3D[:,:,0]) print("check if reflectance grid is build correctly. Reflectance at " + \ "(0.29442218, 0.04258016): " + str(rbs( 0.29442218, 0.04258016 , dx=0)) + "; expected " + str(654.222)) print("check if reflectance grid is build correctly. Reflectance at " + \ "(0.60, 0.1): " + str(rbs( 0.6, 0.1 , dx=0)) + "; expected " + str(818.945)) print("check if reflectance grid is build correctly. Reflectance at " + \ "(0.60, 0.01): " + str(rbs( 0.6, 0.01 , dx=0)) + "; expected " + str(1120.118)) print("check if reflectance grid is build correctly. Reflectance at " + \ "(0.041, 0.0112): " + str(rbsFrom3DGrid(0.041,0.0112, dx=0)) + "; expected " + str(rbs(0.041,0.0112, dx=0))) reflectanceError = ReflectanceError(Vss, BVFs, reflectanceGrid3D) # reflectance at point (0.149, 0.062) exampleReflectance = [ 312.09769118, 327.67866117, 336.20970583, 343.10626114, 344.77202411, 322.68062559, 280.01722363, 265.45441892, 292.59760112, 294.58676191, 256.88475134, 296.10442177, 388.44056814, 422.99479049, 425.6602738 , 426.88353542, 427.09604971, 413.01675144, 424.50008822, 423.78125335, 422.71219033, 421.51283858, 421.32106797] reflectanceError.setReflectanceToMatch(exampleReflectance) print("check if reflectance error value is correctly computed. Expected: 0, actual: " + str(reflectanceError.f([1.149252955326078018e-01, 6.191087469731736126e-02]))) print("check if derivate of reflectance error value is correctly computed. Expected: close to [0,0], actual: " + str(reflectanceError.df([1.149252955326078018e-01, 6.191087469731736126e-02]))) # reflectance at point (0.6, 0.01) exampleReflectance2 = [1.120118468811987896e+03, 1.131318509608477598e+03, 1.137543168078165081e+03, 1.141159082538952589e+03, 1.139712500764461311e+03, 1.120915174518419235e+03, 1.081868496849364647e+03, 1.064766763752600582e+03, 1.084980179526470238e+03, 1.081898683840950525e+03, 1.038063983415024040e+03, 1.073340106470938281e+03, 1.117816020857113244e+03, 1.122940442828165033e+03, 1.116782763468411758e+03, 1.108377361568318520e+03, 1.104745713217586172e+03, 1.098370113333694690e+03, 1.093525893516002043e+03, 1.089713462981120756e+03, 1.085745981406935471e+03, 1.087101893762737973e+03, 1.083097058199290814e+03] reflectanceError.setReflectanceToMatch(exampleReflectance2) print("check if reflectance error value is correctly computed (second one). Expected: 0, actual: " + str(reflectanceError.f([0.6, 0.01]))) print("check if derivate of reflectance error value is correctly computed (second one). Expected: close to [0,0], actual: " + str(reflectanceError.df([0.6, 0.01]))) # check if grid interpolation looks good. #%% grid_x, grid_y = np.mgrid[0.04:0.6:100j, 0.01:0.1:100j] grid_z = rbs.ev(grid_x, grid_y) fig = plt.figure(0) ax = fig.gca(projection='3d') surf = ax.plot_surface(grid_x, grid_y, grid_z, cmap=cm.jet, linewidth=1, antialiased=True) ax.scatter(testParameters[:,1], testParameters[:,0], testReflectances[:,0]) ax.set_zlim3d(np.min(grid_z), np.max(grid_z)) fig.colorbar(surf) plt.show() # plt.imshow(grid_z.T, extent=(0.04,0.6,0.01,0.1))