diff --git a/Modules/Biophotonics/CMakeLists.txt b/Modules/Biophotonics/CMakeLists.txt new file mode 100644 index 0000000000..8b01549c66 --- /dev/null +++ b/Modules/Biophotonics/CMakeLists.txt @@ -0,0 +1,5 @@ +MITK_CREATE_MODULE(mitkBiophotonics +) +if( MITK_USE_Python ) + ADD_SUBDIRECTORY(python) +endif() diff --git a/Modules/Biophotonics/python/TechtreePython.txt b/Modules/Biophotonics/python/TechtreePython.txt new file mode 100644 index 0000000000..c4ba1c784b --- /dev/null +++ b/Modules/Biophotonics/python/TechtreePython.txt @@ -0,0 +1,13 @@ +To execute the python scripts in this folder you will need: + +a working executable to gpumcml (e.g. gpumcml.sm_20) somewhere on your computer. + +Python libraries: +pymiecoated +numpy +scipy 0.15 +sklearn + +Other: +graphviz + diff --git a/Modules/Biophotonics/python/data/colonTemplate.mci b/Modules/Biophotonics/python/data/colonTemplate.mci new file mode 100644 index 0000000000..dac8a3a55d --- /dev/null +++ b/Modules/Biophotonics/python/data/colonTemplate.mci @@ -0,0 +1,24 @@ +#### +# Template of input files for Monte Carlo simulation (mcml). +# Anything in a line after "#" is ignored as comments. +# Space lines are also ignored. +# Lengths are in cm, mua and mus are in 1/cm. +#### + +1.0 # file version +1 # number of runs + +### Specify data for run 1 +outputFilename A # output filename, ASCII/Binary +c_ph # No. of photons +0.01 0.01 # dz, dr +150 100 30 # No. of dz, dr & da. + +4 # No. of layers +# n mua mus g d # One line for each layer +1.0 # n for medium above. +1.0 0 0 1 0.01 # layer 0: air +1.38 muc_c_ua muc_c_us muc_c_g muc_c_d # layer 1: mucosa +1.36 sm_c_ua sm_c_us sm_c_g 0.05 # layer 2: submucosa +1.36 mus_c_ua mus_c_us 0.96 0.09 # layer 3: muscle layer +1.0 # n for medium below. diff --git a/Modules/Biophotonics/python/data/haemoglobin.txt b/Modules/Biophotonics/python/data/haemoglobin.txt new file mode 100644 index 0000000000..b5e87f7035 --- /dev/null +++ b/Modules/Biophotonics/python/data/haemoglobin.txt @@ -0,0 +1,378 @@ +lambda Hb02 Hb +nm cm-1/M cm-1/M +250 106112 112736 +252 105552 112736 +254 107660 112736 +256 109788 113824 +258 112944 115040 +260 116376 116296 +262 120188 117564 +264 124412 118876 +266 128696 120208 +268 133064 121544 +270 136068 122880 +272 137232 123096 +274 138408 121952 +276 137424 120808 +278 135820 119840 +280 131936 118872 +282 127720 117628 +284 122280 114820 +286 116508 112008 +288 108484 107140 +290 104752 98364 +292 98936 91636 +294 88136 85820 +296 79316 77100 +298 70884 69444 +300 65972 64440 +302 63208 61300 +304 61952 58828 +306 62352 56908 +308 62856 57620 +310 63352 59156 +312 65972 62248 +314 69016 65344 +316 72404 68312 +318 75536 71208 +320 78752 74508 +322 82256 78284 +324 85972 82060 +326 89796 85592 +328 93768 88516 +330 97512 90856 +332 100964 93192 +334 103504 95532 +336 104968 99792 +338 106452 104476 +340 107884 108472 +342 109060 110996 +344 110092 113524 +346 109032 116052 +348 107984 118752 +350 106576 122092 +352 105040 125436 +354 103696 128776 +356 101568 132120 +358 97828 133632 +360 94744 134940 +362 92248 136044 +364 89836 136972 +366 88484 137900 +368 87512 138856 +370 88176 139968 +372 91592 141084 +374 95140 142196 +376 98936 143312 +378 103432 144424 +380 109564 145232 +382 116968 145232 +384 125420 148668 +386 135132 153908 +388 148100 159544 +390 167748 167780 +392 189740 180004 +394 212060 191540 +396 231612 202124 +398 248404 212712 +400 266232 223296 +402 284224 236188 +404 308716 253368 +406 354208 270548 +408 422320 287356 +410 466840 303956 +412 500200 321344 +414 524280 342596 +416 521880 363848 +418 515520 385680 +420 480360 407560 +422 431880 429880 +424 376236 461200 +426 326032 481840 +428 283112 500840 +430 246072 528600 +432 214120 552160 +434 165332 552160 +436 132820 547040 +438 119140 501560 +440 102580 413280 +442 92780 363240 +444 81444 282724 +446 76324 237224 +448 67044 173320 +450 62816 103292 +452 58864 62640 +454 53552 36170 +456 49496 30698.8 +458 47496 25886.4 +460 44480 23388.8 +462 41320 20891.2 +464 39807.2 19260.8 +466 37073.2 18142.4 +468 34870.8 17025.6 +470 33209.2 16156.4 +472 31620 15310 +474 30113.6 15048.4 +476 28850.8 14792.8 +478 27718 14657.2 +480 26629.2 14550 +482 25701.6 14881.2 +484 25180.4 15212.4 +486 24669.6 15543.6 +488 24174.8 15898 +490 23684.4 16684 +492 23086.8 17469.6 +494 22457.6 18255.6 +496 21850.4 19041.2 +498 21260 19891.2 +500 20932.8 20862 +502 20596.4 21832.8 +504 20418 22803.6 +506 19946 23774.4 +508 19996 24745.2 +510 20035.2 25773.6 +512 20150.4 26936.8 +514 20429.2 28100 +516 21001.6 29263.2 +518 22509.6 30426.4 +520 24202.4 31589.6 +522 26450.4 32851.2 +524 29269.2 34397.6 +526 32496.4 35944 +528 35990 37490 +530 39956.8 39036.4 +532 43876 40584 +534 46924 42088 +536 49752 43592 +538 51712 45092 +540 53236 46592 +542 53292 48148 +544 52096 49708 +546 49868 51268 +548 46660 52496 +550 43016 53412 +552 39675.2 54080 +554 36815.2 54520 +556 34476.8 54540 +558 33456 54164 +560 32613.2 53788 +562 32620 52276 +564 33915.6 50572 +566 36495.2 48828 +568 40172 46948 +570 44496 45072 +572 49172 43340 +574 53308 41716 +576 55540 40092 +578 54728 38467.6 +580 50104 37020 +582 43304 35676.4 +584 34639.6 34332.8 +586 26600.4 32851.6 +588 19763.2 31075.2 +590 14400.8 28324.4 +592 10468.4 25470 +594 7678.8 22574.8 +596 5683.6 19800 +598 4504.4 17058.4 +600 3200 14677.2 +602 2664 13622.4 +604 2128 12567.6 +606 1789.2 11513.2 +608 1647.6 10477.6 +610 1506 9443.6 +612 1364.4 8591.2 +614 1222.8 7762 +616 1110 7344.8 +618 1026 6927.2 +620 942 6509.6 +622 858 6193.2 +624 774 5906.8 +626 707.6 5620 +628 658.8 5366.8 +630 610 5148.8 +632 561.2 4930.8 +634 512.4 4730.8 +636 478.8 4602.4 +638 460.4 4473.6 +640 442 4345.2 +642 423.6 4216.8 +644 405.2 4088.4 +646 390.4 3965.08 +648 379.2 3857.6 +650 368 3750.12 +652 356.8 3642.64 +654 345.6 3535.16 +656 335.2 3427.68 +658 325.6 3320.2 +660 319.6 3226.56 +662 314 3140.28 +664 308.4 3053.96 +666 302.8 2967.68 +668 298 2881.4 +670 294 2795.12 +672 290 2708.84 +674 285.6 2627.64 +676 282 2554.4 +678 279.2 2481.16 +680 277.6 2407.92 +682 276 2334.68 +684 274.4 2261.48 +686 272.8 2188.24 +688 274.4 2115 +690 276 2051.96 +692 277.6 2000.48 +694 279.2 1949.04 +696 282 1897.56 +698 286 1846.08 +700 290 1794.28 +702 294 1741 +704 298 1687.76 +706 302.8 1634.48 +708 308.4 1583.52 +710 314 1540.48 +712 319.6 1497.4 +714 325.2 1454.36 +716 332 1411.32 +718 340 1368.28 +720 348 1325.88 +722 356 1285.16 +724 364 1244.44 +726 372.4 1203.68 +728 381.2 1152.8 +730 390 1102.2 +732 398.8 1102.2 +734 407.6 1102.2 +736 418.8 1101.76 +738 432.4 1100.48 +740 446 1115.88 +742 459.6 1161.64 +744 473.2 1207.4 +746 487.6 1266.04 +748 502.8 1333.24 +750 518 1405.24 +752 533.2 1515.32 +754 548.4 1541.76 +756 562 1560.48 +758 574 1560.48 +760 586 1548.52 +762 598 1508.44 +764 610 1459.56 +766 622.8 1410.52 +768 636.4 1361.32 +770 650 1311.88 +772 663.6 1262.44 +774 677.2 1213 +776 689.2 1163.56 +778 699.6 1114.8 +780 710 1075.44 +782 720.4 1036.08 +784 730.8 996.72 +786 740 957.36 +788 748 921.8 +790 756 890.8 +792 764 859.8 +794 772 828.8 +796 786.4 802.96 +798 807.2 782.36 +800 816 761.72 +802 828 743.84 +804 836 737.08 +806 844 730.28 +808 856 723.52 +810 864 717.08 +812 872 711.84 +814 880 706.6 +816 887.2 701.32 +818 901.6 696.08 +820 916 693.76 +822 930.4 693.6 +824 944.8 693.48 +826 956.4 693.32 +828 965.2 693.2 +830 974 693.04 +832 982.8 692.92 +834 991.6 692.76 +836 1001.2 692.64 +838 1011.6 692.48 +840 1022 692.36 +842 1032.4 692.2 +844 1042.8 691.96 +846 1050 691.76 +848 1054 691.52 +850 1058 691.32 +852 1062 691.08 +854 1066 690.88 +856 1072.8 690.64 +858 1082.4 692.44 +860 1092 694.32 +862 1101.6 696.2 +864 1111.2 698.04 +866 1118.4 699.92 +868 1123.2 701.8 +870 1128 705.84 +872 1132.8 709.96 +874 1137.6 714.08 +876 1142.8 718.2 +878 1148.4 722.32 +880 1154 726.44 +882 1159.6 729.84 +884 1165.2 733.2 +886 1170 736.6 +888 1174 739.96 +890 1178 743.6 +892 1182 747.24 +894 1186 750.88 +896 1190 754.52 +898 1194 758.16 +900 1198 761.84 +902 1202 765.04 +904 1206 767.44 +906 1209.2 769.8 +908 1211.6 772.16 +910 1214 774.56 +912 1216.4 776.92 +914 1218.8 778.4 +916 1220.8 778.04 +918 1222.4 777.72 +920 1224 777.36 +922 1225.6 777.04 +924 1227.2 776.64 +926 1226.8 772.36 +928 1224.4 768.08 +930 1222 763.84 +932 1219.6 752.28 +934 1217.2 737.56 +936 1215.6 722.88 +938 1214.8 708.16 +940 1214 693.44 +942 1213.2 678.72 +944 1212.4 660.52 +946 1210.4 641.08 +948 1207.2 621.64 +950 1204 602.24 +952 1200.8 583.4 +954 1197.6 568.92 +956 1194 554.48 +958 1190 540.04 +960 1186 525.56 +962 1182 511.12 +964 1178 495.36 +966 1173.2 473.32 +968 1167.6 451.32 +970 1162 429.32 +972 1156.4 415.28 +974 1150.8 402.28 +976 1144 389.288 +978 1136 374.944 +980 1128 359.656 +982 1120 344.372 +984 1112 329.084 +986 1102.4 313.796 +988 1091.2 298.508 +990 1080 283.22 +992 1068.8 267.932 +994 1057.6 252.648 +996 1046.4 237.36 +998 1035.2 222.072 +1000 1024 206.784 \ No newline at end of file diff --git a/Modules/Biophotonics/python/generateNoisyRandomSpectra.py b/Modules/Biophotonics/python/generateNoisyRandomSpectra.py new file mode 100644 index 0000000000..8f48790d92 --- /dev/null +++ b/Modules/Biophotonics/python/generateNoisyRandomSpectra.py @@ -0,0 +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 = '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, 8)) + +print('start simulations...') + +#%% start program logic +start = time.time() + + + +for i in range(nrSimulations): + + 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 j, wavelength in enumerate(wavelengths): + + reflectanceValue = mch.runOneSimulation( + wavelength, eHbO2, eHb, + 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 + + + + +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 new file mode 100644 index 0000000000..98b5382421 --- /dev/null +++ b/Modules/Biophotonics/python/generateRandomSpectra.py @@ -0,0 +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): + + 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 j, wavelength in enumerate(wavelengths): + + + reflectanceValue = mch.runOneSimulation( + wavelength, eHbO2, eHb, + 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 + + + + +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/generateReflectanceSpectra.py b/Modules/Biophotonics/python/generateReflectanceSpectra.py new file mode 100644 index 0000000000..4ed4bbd364 --- /dev/null +++ b/Modules/Biophotonics/python/generateReflectanceSpectra.py @@ -0,0 +1,87 @@ +import time +import itertools +import datetime +import os + +import numpy as np + +import matplotlib.pyplot as plt + +import monteCarloHelper as mch + +import setupSimulation as setup + +""" +script generating: +_____ +1. + All the monte carlo simulations +2. + Reflectance over wavelengths for the given parameter ranges (assuming a + infinitly wide, uniform light source) + +TODO: +_____ + +""" + +# 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() + +reflectances = np.zeros((nrSamples, len(wavelengths))) + +print('start simulations...') + + +#%% start program logic +start = time.time() + +currentSimulation = 0 + +params = itertools.product(BVFs, Vss, ds, SaO2s, rs) +paramsList = list(itertools.product(BVFs, Vss, ds, SaO2s, rs)) + +for BVF, Vs, d, SaO2, r in params: + j = 0 + print('starting simulation ' + str(currentSimulation) + ' of ' + str(nrSamples)) + for wavelength in wavelengths: + + + reflectanceValue = mch.runOneSimulation( + wavelength, eHbO2, eHb, + 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[currentSimulation, j] = reflectanceValue + j +=1 + + currentSimulation += 1 + + + + +infile.close() + +# save the reflectance results! +now = datetime.datetime.now().strftime("%Y%B%d%I:%M%p") +np.save(outfolderRS + now + "reflectances" + str(photons) + "photons", reflectances) +np.save(outfolderRS + now + "paramters", paramsList) + +end = time.time() +print "total time for script: " + str((end - start)) \ No newline at end of file diff --git a/Modules/Biophotonics/python/mieMonteCarlo.py b/Modules/Biophotonics/python/mieMonteCarlo.py new file mode 100644 index 0000000000..3c5c526d03 --- /dev/null +++ b/Modules/Biophotonics/python/mieMonteCarlo.py @@ -0,0 +1,72 @@ +# -*- coding: utf-8 -*- + + +import math +from pymiecoated import Mie + +#%% setup and create derived paramters + +def mieMonteCarlo(wavelength = 450*10**-9, r = 300*10**-9, Vs = 0.1, nParticle = 1.46, nMedium = 1.36): + """ + Calculate the scattering parameters relevant for monte carlo simulation. + + Needs pymiecoated: https://code.google.com/p/pymiecoated/ + These are calculate scattering coefficient [1/cm] and anisotropy factor for given: + + Args + ____ + wavelength: + wavelength of the incident light [m] + r: + radius of the particle [m] + Vs: + volume fraction of scattering particles + nParticle: + refractive index of the particle that the light wave is scattered on + (default value is the refractive index of collagen) + nMedium: + refractive index of the surronding medium (default is that of colonic + mucosal tissue) + + Returns: + ____ + {'us', 'g'}: + scattering coefficient us [1/m] and anisotropy factor g + + TODO: + _____ + Additional input parameter specifying a FWHM for the wavelength to simulate the + scattering for a broad filter + """ + + # create derived parameters + sizeParamter = 2 * math.pi * r / wavelength + nRelative = nParticle / nMedium + + #%% execute mie and create derived parameters + + mie = Mie(x=sizeParamter, m=complex(nRelative,0.0)) # 0.0 complex for no attenuation + + A = math.pi * r**2 # geometrical cross sectional area + cs = mie.qsca() * A # scattering cross section + + us = Vs / (4/3 * r**3 * math.pi) * cs # scattering coefficient [m⁻1] + + return {'us': us, 'g': mie.asy()} + +def mieMonteCarlo_FWHM(wavelength = 450*10**-9, r = 300*10**-9, Vs = 0.1, nParticle = 1.46, nMedium = 1.36, FWHM = 0): + vl = mieMonteCarlo(wavelength - FWHM / 2, r, Vs, nParticle, nMedium) + vr = mieMonteCarlo(wavelength + FWHM / 2, r, Vs, nParticle, nMedium) + vc = mieMonteCarlo(wavelength, r, Vs, nParticle, nMedium) + + us = (0.5 * vl['us'] + vc['us'] + 0.5 * vr['us']) / 2 + g = (0.5 * vl['g'] + vc['g'] + 0.5 * vr['g']) / 2 + return {'us': us, 'g': g} + +if __name__ == "__main__": + + usg = mieMonteCarlo(wavelength = 450*10**-9, r = 300*10**-9, Vs = 0.1, nParticle = 1.46, nMedium = 1.36) + + #%% print results + print("Scattering coefficient in [1/cm]: " + str(usg['us'] * 10**-2) + "; expected: 589.36") + print("Anisotropy factor: " + str(usg['g']) + "; expected: 0.88") diff --git a/Modules/Biophotonics/python/monteCarloHelper.py b/Modules/Biophotonics/python/monteCarloHelper.py new file mode 100644 index 0000000000..1e0975d74e --- /dev/null +++ b/Modules/Biophotonics/python/monteCarloHelper.py @@ -0,0 +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, 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 = ("./" + 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 new file mode 100644 index 0000000000..6a11a8c649 --- /dev/null +++ b/Modules/Biophotonics/python/optimization.py @@ -0,0 +1,154 @@ +# -*- coding: utf-8 -*- +""" +Created on Sat Feb 07 23:42:41 2015 + +@author: Seb +""" + +import numpy as np +import matplotlib.pyplot as plt +from scipy.interpolate import RectBivariateSpline +from scipy.optimize import minimize +from reflectanceError import ReflectanceError +from mpl_toolkits.mplot3d import Axes3D +from matplotlib import cm + +import setupData + + +#%% load data + +dataFolder = "data/output/" + +trainingParameters, trainingReflectances, testParameters, testReflectances = \ + setupData.setupTwoDimensionalData(dataFolder) + + +BVFs = np.unique(trainingParameters[:,0]) +Vss = np.unique(trainingParameters[:,1]) + + +#%% build optimization function from rectangular reflectance grid + +reflectanceGrid3D = np.reshape(trainingReflectances, (len(BVFs), len(Vss), trainingReflectances.shape[1])) + +functionToMinimize = ReflectanceError(BVFs, Vss, reflectanceGrid3D) + + +#%% do optimization + +absErrors = np.zeros_like(testParameters) + + +for idx, (testParameter, testReflectance) in enumerate(zip(testParameters, testReflectances)): + functionToMinimize.setReflectanceToMatch(testReflectance) + minimization = minimize(functionToMinimize.f, [0.05, 0.3], method="Nelder-Mead") + # interpolation extrapolates with constant values. Since the minimization function is + # covex, we can just crop it to the bounds + estimatedBVF= np.clip(minimization.x[0], min(BVFs), max(BVFs)) + estimatedVs = np.clip(minimization.x[1], min(Vss), max(Vss)) + + absErrors[idx,:] = np.abs([estimatedBVF, estimatedVs] - testParameter) + + +#%% test + +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))) + + + +#%% test routines, put these in a unit test asap! +if __name__ == "__main__": + + rbs = RectBivariateSpline(BVFs, + Vss, + np.reshape(trainingReflectances[:,0],(100,100))) + + + + rbsFrom3DGrid = RectBivariateSpline(BVFs, + Vss, + reflectanceGrid3D[:,:,0]) + + print("check if reflectance grid is build correctly. Reflectance at " + \ + "(0.04258016, 0.29442218): " + str(rbs( 0.04258016 , 0.29442218, dx=0)) + "; expected " + str(654.222)) + + print("check if reflectance grid is build correctly. Reflectance at " + \ + "(0.1, 0.60): " + str(rbs( 0.1 , 0.6, dx=0)) + "; expected " + str(818.945)) + + print("check if reflectance grid is build correctly. Reflectance at " + \ + "(0.01, 0.60): " + str(rbs(0.01 , 0.6, dx=0)) + "; expected " + str(1120.118)) + + + print("check if reflectance grid is build correctly. Reflectance at " + \ + "(0.0112, 0.041): " + str(rbsFrom3DGrid(0.0112, 0.041, dx=0)) + "; expected " + str(rbs(0.0112, 0.041, dx=0))) + + reflectanceError = ReflectanceError(BVFs, Vss, reflectanceGrid3D) + + # reflectance at point (0.062, 0.149) + 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([6.191087469731736126e-02, 1.149252955326078018e-01]))) + + print("check if derivate of reflectance error value is correctly computed. Expected: close to [0,0], actual: " + + str(reflectanceError.df([6.191087469731736126e-02, 1.149252955326078018e-01]))) + + # reflectance at point (0.01, 0.6) + 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.01, 0.6]))) + + print("check if derivate of reflectance error value is correctly computed (second one). Expected: close to [0,0], actual: " + + str(reflectanceError.df([0.01, 0.6]))) + + + functionToMinimize.setReflectanceToMatch(trainingReflectances[0,:]) + absError = minimize(functionToMinimize.f, [0.05, 0.3], method="Nelder-Mead") + #absError = minimize(functionToMinimize.f, [0.05, 0.3], method="BFGS", jac=functionToMinimize.df) + + print(" function evaluation yields: " + str(absError.x) + " (unclipped) with absolute error: " + str(functionToMinimize.f(absError.x))) + print(" expected evaluation result: " + str([0.01, 0.04]) + " real result error: " + str(functionToMinimize.f([0.01, 0.04]))) + + + + + # check if grid interpolation looks good. + #%% + grid_x, grid_y = np.mgrid[min(BVFs):max(BVFs):100j, min(Vss):max(Vss):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[:,0], testParameters[:,1], 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)) + diff --git a/Modules/Biophotonics/python/randomForest.py b/Modules/Biophotonics/python/randomForest.py new file mode 100644 index 0000000000..5e84957ab3 --- /dev/null +++ b/Modules/Biophotonics/python/randomForest.py @@ -0,0 +1,52 @@ +# -*- coding: utf-8 -*- +""" +Created on Fri Feb 6 10:49:45 2015 + +@author: wirkert +""" + + +import numpy as np +import matplotlib.pyplot as plt +from sklearn import tree +from sklearn.ensemble import RandomForestRegressor + +import setupData + + + +# additional things that could be checked in this study: +# 1. band selection results +# 2. effect of normalizations +# 3. parameter study +# 4. optimal image quotient + +#%% initialize + + +# the folder with the reflectance spectra +dataFolder = 'data/output/' + +# load data +trainingParameters, trainingReflectances, testParameters, testReflectances = \ + setupData.setupTwoDimensionalData(dataFolder) + + +#%% train forest + +rf = RandomForestRegressor(50, max_depth=8, max_features="log2", n_jobs=5) +rf.fit(trainingReflectances, trainingParameters) + +#with open("iris.dot", 'w') as f: +# f = tree.export_graphviz(rf, out_file=f) + +#%% test + + +# predict test reflectances and get absolute errors. +absErrors = np.abs(rf.predict(testReflectances) - testParameters) + +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))) diff --git a/Modules/Biophotonics/python/reflectanceError.py b/Modules/Biophotonics/python/reflectanceError.py new file mode 100644 index 0000000000..fcf4b63926 --- /dev/null +++ b/Modules/Biophotonics/python/reflectanceError.py @@ -0,0 +1,75 @@ +# -*- coding: utf-8 -*- +""" +Created on Tue Feb 10 15:50:13 2015 + +@author: wirkert +""" + +from scipy.interpolate import RectBivariateSpline +import numpy as np + +class ReflectanceError: + """ + calculate the fit error of a given reflectance and examples of reflectances. + """ + + def __init__(self, BVFs, Vss, reflectances): + """ + intialize the reflectance error calculation. + + Arguments + _________ + Vss: + array storing all the different values Vs (Volume fraction of scattering particles) + BVFs: + array storing all the differnet values for BVF (Blood Volume Fraction) + reflectances: + grid of dimension Vss.shape[0] x BVFs.shape[0] x wavelenghts with reflectance values + """ + nr_wavelengths = reflectances.shape[2] + + # create array of bivariate splines, one for each wavelength + self.RBSs = np.empty(nr_wavelengths, dtype=object) + + for i in np.arange(nr_wavelengths): + self.RBSs[i] = RectBivariateSpline(BVFs, + Vss, + reflectances[:,:,i]) + + def setReflectanceToMatch(self, reflectance): + """ + set the baseline reflectance. the l2 distance to this reflectance will be output + """ + self.reflectance = reflectance + + def reflectanceForX(self, xy): + result = np.zeros(len(self.reflectance)) + for idx, rbsi in enumerate(self.RBSs): + result[idx] = rbsi(xy[0], xy[1]) + return result + + def f(self, xy): + """ + return l2 norm of reflectance implied by x and self.reflectances + """ + quadraticError = 0 + # build quadratic error of differences between reflectances and + for refli, rbsi in zip(self.reflectance, self.RBSs): + quadraticError += (refli - rbsi(xy[0], xy[1], grid = False))**2 + + return quadraticError + + def df(self, xy): + """ + return derivate of l2 norm of reflectance implied by x and self.reflectances + """ + quadraticErrorDerivate = [0,0] + for refli, rbsi in zip(self.reflectance, self.RBSs): + u = 2*(rbsi(xy[0], xy[1], grid = False) - refli) + # dx + quadraticErrorDerivate[0] += \ + u*rbsi(xy[0], xy[1], dx=1, grid = False) + #dy + quadraticErrorDerivate[1] += \ + u*rbsi(xy[0], xy[1], dy=1, grid = False) + return quadraticErrorDerivate diff --git a/Modules/Biophotonics/python/setupData.py b/Modules/Biophotonics/python/setupData.py new file mode 100644 index 0000000000..b61517bf49 --- /dev/null +++ b/Modules/Biophotonics/python/setupData.py @@ -0,0 +1,27 @@ +# -*- coding: utf-8 -*- +""" +Created on Fri Feb 13 16:13:11 2015 + +functions storing training / testdata combinations to ensure consistent usage +of datasets (same combinations / normalizations / ...) + +@author: wirkert +""" + +import numpy as np +import monteCarloHelper as mch + +def setupTwoDimensionalData(dataFolder): + trainingParameters = np.load(dataFolder + "2015February0511:02PMparamters2D.npy") + trainingReflectances = np.load(dataFolder + "2015February0511:02PMreflectances2D.npy") + # trainingReflectances = mch.normalizeImageQuotient(trainingReflectances, iqBand=4) + + testParameters = np.load(dataFolder + "2015February1107:43PMparamtersRandom2D.npy") + testReflectances = np.load(dataFolder + "2015February1107:43PMreflectancesRandom2D.npy") + # testReflectances = mch.normalizeImageQuotient(testReflectances, iqBand=4) + + return trainingParameters, trainingReflectances, testParameters, testReflectances + + + + diff --git a/Modules/Biophotonics/python/setupSimulation.py b/Modules/Biophotonics/python/setupSimulation.py new file mode 100644 index 0000000000..874e898f71 --- /dev/null +++ b/Modules/Biophotonics/python/setupSimulation.py @@ -0,0 +1,75 @@ +# -*- coding: utf-8 -*- +""" +Created on Fri Feb 6 15:32:26 2015 + +@author: wirkert +""" + + + +import numpy as np +from scipy.interpolate import interp1d +from scipy.stats import norm + +def setupNormalSimulation(): + #%% intialization + print('initialize...') + + # evaluated parameters for MC + + # blood volume fraction range + samplesBVF = 20 + BVFs = np.linspace(0.01, 0.20, samplesBVF) + # scattering particles volume fraction range + samplesVs = 20 + Vss = np.linspace(0.04, 0.60, samplesVs) + # thickness of mucosa [m] range + samplesD = 5 + ds = np.linspace(250 * 10**-6, 735 * 10**-6, samplesD) + #ds = np.linspace(500 * 10**-6, 500 * 10**-6, samplesD) + # haemoglobin saturation + samplesSaO2 = 6 + SaO2s= np.linspace(0.3, 0.9, samplesSaO2) + # radius of scattering particles [m] + samplesR = 1 + rs = np.linspace(0.4 * 10**-6, 0.4 * 10**-6, samplesR) + # number of photons for one MC simulation run + photons = 1 * 10**6 + # the wavelengths [m] for which the reflectance spectrum shall be evaluated + #wavelengths = np.linspace(470,700,23) * 10**-9 + wavelengths = np.array([470, 480, 511, 560, 580, 600, 660, 700]) * 10**-9 + + # The full width at half maximum [m] of the used imaging systems filters + FWHM = 20 * 10**-9 + + print('create reference data...') + #reference data + + # table with wavelength at 1st row, + # HbO2 molar extinction coefficient [cm**-1/(moles/l)] at 2nd row, + # Hb molar extinction coefficient [cm**-1/(moles/l)] at 3rd row + haemoLUT = np.loadtxt("data/haemoglobin.txt", skiprows=2) + # we calculate everything in [m] instead of [nm] and [1/cm] + haemoLUT[:,0] = haemoLUT[:,0] * 10**-9 + haemoLUT[:,1:2] = haemoLUT[:,1:2] * 10**2 + eHbO2 = interp1d(haemoLUT[:,0], haemoLUT[:,1]) + eHb = interp1d(haemoLUT[:,0], haemoLUT[:,2]) + # to account for the FWHM of the used filters, compute convolution + # see http://en.wikipedia.org/wiki/Full_width_at_half_maximum + filterResponse = norm(loc = 0, scale = FWHM / 2.355) + # parse at 20 locations + x = np.linspace(filterResponse.ppf(0.01), + filterResponse.ppf(0.99), 20) + filterResponse_table = filterResponse.pdf(x) + # TODO verify if this normalization is correct! + # filterResponse_table = filterResponse_table / sum(filterResponse_table) + haemoLUT[:, 1] = np.convolve(haemoLUT[:, 1], filterResponse_table, 'same') + haemoLUT[:, 2] = np.convolve(haemoLUT[:, 2], filterResponse_table, 'same') + #plt.plot(haemoLUT[:, 0], haemoLUT[:, 1], '--', linewidth=2) + #plt.show() + + + nrSamples = samplesBVF * samplesD * samplesR * samplesSaO2 * samplesVs + + + return BVFs, Vss, ds, SaO2s, rs, nrSamples, photons, wavelengths, FWHM, eHbO2, eHb