diff options
-rw-r--r-- | chromaticity.py | 18 | ||||
-rw-r--r-- | ciexyz64_1.csv | 471 | ||||
-rw-r--r-- | geometry.py | 48 | ||||
-rw-r--r-- | materials.py | 10 | ||||
-rw-r--r-- | models/Colbert_HighRes_Brow.STL | bin | 0 -> 13910384 bytes | |||
-rw-r--r-- | models/Colbert_HighRes_Smile.STL | bin | 0 -> 8699284 bytes | |||
-rw-r--r-- | models/film.stl | 86 | ||||
-rw-r--r-- | models/liberty.stl | bin | 0 -> 853084 bytes | |||
-rw-r--r-- | photon_map.py | 86 | ||||
-rw-r--r-- | sample.py (renamed from photon.py) | 0 | ||||
-rw-r--r-- | solid.py | 19 | ||||
-rw-r--r-- | src/intersect.h | 6 | ||||
-rw-r--r-- | src/kernel.cu | 268 | ||||
-rw-r--r-- | src/materials.h | 45 | ||||
-rw-r--r-- | src/physical_constants.h | 7 | ||||
-rw-r--r-- | src/random.h | 22 |
16 files changed, 1009 insertions, 77 deletions
diff --git a/chromaticity.py b/chromaticity.py new file mode 100644 index 0000000..18ec28f --- /dev/null +++ b/chromaticity.py @@ -0,0 +1,18 @@ +import numpy as np + +f = open('ciexyz64_1.csv') + +color_map = [] +for line in f: + color_map.append([float(s) for s in line.split(',')]) + +f.close() + +color_map = np.array(color_map) + +def map_wavelength(wavelength): + r = np.interp(wavelength, color_map[:,0], color_map[:,1]) + g = np.interp(wavelength, color_map[:,0], color_map[:,2]) + b = np.interp(wavelength, color_map[:,0], color_map[:,3]) + + return r, g, b diff --git a/ciexyz64_1.csv b/ciexyz64_1.csv new file mode 100644 index 0000000..f437db7 --- /dev/null +++ b/ciexyz64_1.csv @@ -0,0 +1,471 @@ +360,0.000000122200,0.000000013398,0.000000535027
+361,0.000000185138,0.000000020294,0.000000810720
+362,0.000000278830,0.000000030560,0.000001221200
+363,0.000000417470,0.000000045740,0.000001828700
+364,0.000000621330,0.000000068050,0.000002722200
+365,0.000000919270,0.000000100650,0.000004028300
+366,0.000001351980,0.000000147980,0.000005925700
+367,0.000001976540,0.000000216270,0.000008665100
+368,0.000002872500,0.000000314200,0.000012596000
+369,0.000004149500,0.000000453700,0.000018201000
+370,0.000005958600,0.000000651100,0.000026143700
+371,0.000008505600,0.000000928800,0.000037330000
+372,0.000012068600,0.000001317500,0.000052987000
+373,0.000017022600,0.000001857200,0.000074764000
+374,0.000023868000,0.000002602000,0.000104870000
+375,0.000033266000,0.000003625000,0.000146220000
+376,0.000046087000,0.000005019000,0.000202660000
+377,0.000063472000,0.000006907000,0.000279230000
+378,0.000086892000,0.000009449000,0.000382450000
+379,0.000118246000,0.000012848000,0.000520720000
+380,0.000159952000,0.000017364000,0.000704776000
+381,0.000215080000,0.000023327000,0.000948230000
+382,0.000287490000,0.000031150000,0.001268200000
+383,0.000381990000,0.000041350000,0.001686100000
+384,0.000504550000,0.000054560000,0.002228500000
+385,0.000662440000,0.000071560000,0.002927800000
+386,0.000864500000,0.000093300000,0.003823700000
+387,0.001121500000,0.000120870000,0.004964200000
+388,0.001446160000,0.000155640000,0.006406700000
+389,0.001853590000,0.000199200000,0.008219300000
+390,0.002361600000,0.000253400000,0.010482200000
+391,0.002990600000,0.000320200000,0.013289000000
+392,0.003764500000,0.000402400000,0.016747000000
+393,0.004710200000,0.000502300000,0.020980000000
+394,0.005858100000,0.000623200000,0.026127000000
+395,0.007242300000,0.000768500000,0.032344000000
+396,0.008899600000,0.000941700000,0.039802000000
+397,0.010870900000,0.001147800000,0.048691000000
+398,0.013198900000,0.001390300000,0.059210000000
+399,0.015929200000,0.001674000000,0.071576000000
+400,0.019109700000,0.002004400000,0.086010900000
+401,0.022788000000,0.002386000000,0.102740000000
+402,0.027011000000,0.002822000000,0.122000000000
+403,0.031829000000,0.003319000000,0.144020000000
+404,0.037278000000,0.003880000000,0.168990000000
+405,0.043400000000,0.004509000000,0.197120000000
+406,0.050223000000,0.005209000000,0.228570000000
+407,0.057764000000,0.005985000000,0.263470000000
+408,0.066038000000,0.006833000000,0.301900000000
+409,0.075033000000,0.007757000000,0.343870000000
+410,0.084736000000,0.008756000000,0.389366000000
+411,0.095041000000,0.009816000000,0.437970000000
+412,0.105836000000,0.010918000000,0.489220000000
+413,0.117066000000,0.012058000000,0.542900000000
+414,0.128682000000,0.013237000000,0.598810000000
+415,0.140638000000,0.014456000000,0.656760000000
+416,0.152893000000,0.015717000000,0.716580000000
+417,0.165416000000,0.017025000000,0.778120000000
+418,0.178191000000,0.018399000000,0.841310000000
+419,0.191214000000,0.019848000000,0.906110000000
+420,0.204492000000,0.021391000000,0.972542000000
+421,0.217650000000,0.022992000000,1.038900000000
+422,0.230267000000,0.024598000000,1.103100000000
+423,0.242311000000,0.026213000000,1.165100000000
+424,0.253793000000,0.027841000000,1.224900000000
+425,0.264737000000,0.029497000000,1.282500000000
+426,0.275195000000,0.031195000000,1.338200000000
+427,0.285301000000,0.032927000000,1.392600000000
+428,0.295143000000,0.034738000000,1.446100000000
+429,0.304869000000,0.036654000000,1.499400000000
+430,0.314679000000,0.038676000000,1.553480000000
+431,0.324355000000,0.040792000000,1.607200000000
+432,0.333570000000,0.042946000000,1.658900000000
+433,0.342243000000,0.045114000000,1.708200000000
+434,0.350312000000,0.047333000000,1.754800000000
+435,0.357719000000,0.049602000000,1.798500000000
+436,0.364482000000,0.051934000000,1.839200000000
+437,0.370493000000,0.054337000000,1.876600000000
+438,0.375727000000,0.056822000000,1.910500000000
+439,0.380158000000,0.059399000000,1.940800000000
+440,0.383734000000,0.062077000000,1.967280000000
+441,0.386327000000,0.064737000000,1.989100000000
+442,0.387858000000,0.067285000000,2.005700000000
+443,0.388396000000,0.069764000000,2.017400000000
+444,0.387978000000,0.072218000000,2.024400000000
+445,0.386726000000,0.074704000000,2.027300000000
+446,0.384696000000,0.077272000000,2.026400000000
+447,0.382006000000,0.079979000000,2.022300000000
+448,0.378709000000,0.082874000000,2.015300000000
+449,0.374915000000,0.086000000000,2.006000000000
+450,0.370702000000,0.089456000000,1.994800000000
+451,0.366089000000,0.092947000000,1.981400000000
+452,0.361045000000,0.096275000000,1.965300000000
+453,0.355518000000,0.099535000000,1.946400000000
+454,0.349486000000,0.102829000000,1.924800000000
+455,0.342957000000,0.106256000000,1.900700000000
+456,0.335893000000,0.109901000000,1.874100000000
+457,0.328284000000,0.113835000000,1.845100000000
+458,0.320150000000,0.118167000000,1.813900000000
+459,0.311475000000,0.122932000000,1.780600000000
+460,0.302273000000,0.128201000000,1.745370000000
+461,0.292858000000,0.133457000000,1.709100000000
+462,0.283502000000,0.138323000000,1.672300000000
+463,0.274044000000,0.143042000000,1.634700000000
+464,0.264263000000,0.147787000000,1.595600000000
+465,0.254085000000,0.152761000000,1.554900000000
+466,0.243392000000,0.158102000000,1.512200000000
+467,0.232187000000,0.163941000000,1.467300000000
+468,0.220488000000,0.170362000000,1.419900000000
+469,0.208198000000,0.177425000000,1.370000000000
+470,0.195618000000,0.185190000000,1.317560000000
+471,0.183034000000,0.193025000000,1.262400000000
+472,0.170222000000,0.200313000000,1.205000000000
+473,0.157348000000,0.207156000000,1.146600000000
+474,0.144650000000,0.213644000000,1.088000000000
+475,0.132349000000,0.219940000000,1.030200000000
+476,0.120584000000,0.226170000000,0.973830000000
+477,0.109456000000,0.232467000000,0.919430000000
+478,0.099042000000,0.239025000000,0.867460000000
+479,0.089388000000,0.245997000000,0.818280000000
+480,0.080507000000,0.253589000000,0.772125000000
+481,0.072034000000,0.261876000000,0.728290000000
+482,0.063710000000,0.270643000000,0.686040000000
+483,0.055694000000,0.279645000000,0.645530000000
+484,0.048117000000,0.288694000000,0.606850000000
+485,0.041072000000,0.297665000000,0.570060000000
+486,0.034642000000,0.306469000000,0.535220000000
+487,0.028896000000,0.315035000000,0.502340000000
+488,0.023876000000,0.323335000000,0.471400000000
+489,0.019628000000,0.331366000000,0.442390000000
+490,0.016172000000,0.339133000000,0.415254000000
+491,0.013300000000,0.347860000000,0.390024000000
+492,0.010759000000,0.358326000000,0.366399000000
+493,0.008542000000,0.370001000000,0.344015000000
+494,0.006661000000,0.382464000000,0.322689000000
+495,0.005132000000,0.395379000000,0.302356000000
+496,0.003982000000,0.408482000000,0.283036000000
+497,0.003239000000,0.421588000000,0.264816000000
+498,0.002934000000,0.434619000000,0.247848000000
+499,0.003114000000,0.447601000000,0.232318000000
+500,0.003816000000,0.460777000000,0.218502000000
+501,0.005095000000,0.474340000000,0.205851000000
+502,0.006936000000,0.488200000000,0.193596000000
+503,0.009299000000,0.502340000000,0.181736000000
+504,0.012147000000,0.516740000000,0.170281000000
+505,0.015444000000,0.531360000000,0.159249000000
+506,0.019156000000,0.546190000000,0.148673000000
+507,0.023250000000,0.561180000000,0.138609000000
+508,0.027690000000,0.576290000000,0.129096000000
+509,0.032444000000,0.591500000000,0.120215000000
+510,0.037465000000,0.606741000000,0.112044000000
+511,0.042956000000,0.622150000000,0.104710000000
+512,0.049114000000,0.637830000000,0.098196000000
+513,0.055920000000,0.653710000000,0.092361000000
+514,0.063349000000,0.669680000000,0.087088000000
+515,0.071358000000,0.685660000000,0.082248000000
+516,0.079901000000,0.701550000000,0.077744000000
+517,0.088909000000,0.717230000000,0.073456000000
+518,0.098293000000,0.732570000000,0.069268000000
+519,0.107949000000,0.747460000000,0.065060000000
+520,0.117749000000,0.761757000000,0.060709000000
+521,0.127839000000,0.775340000000,0.056457000000
+522,0.138450000000,0.788220000000,0.052609000000
+523,0.149516000000,0.800460000000,0.049122000000
+524,0.161041000000,0.812140000000,0.045954000000
+525,0.172953000000,0.823330000000,0.043050000000
+526,0.185209000000,0.834120000000,0.040368000000
+527,0.197755000000,0.844600000000,0.037839000000
+528,0.210538000000,0.854870000000,0.035384000000
+529,0.223460000000,0.865040000000,0.032949000000
+530,0.236491000000,0.875211000000,0.030451000000
+531,0.249633000000,0.885370000000,0.028029000000
+532,0.262972000000,0.895370000000,0.025862000000
+533,0.276515000000,0.905150000000,0.023920000000
+534,0.290269000000,0.914650000000,0.022174000000
+535,0.304213000000,0.923810000000,0.020584000000
+536,0.318361000000,0.932550000000,0.019127000000
+537,0.332705000000,0.940810000000,0.017740000000
+538,0.347232000000,0.948520000000,0.016403000000
+539,0.361926000000,0.955600000000,0.015064000000
+540,0.376772000000,0.961988000000,0.013676000000
+541,0.391683000000,0.967540000000,0.012308000000
+542,0.406594000000,0.972230000000,0.011056000000
+543,0.421539000000,0.976170000000,0.009915000000
+544,0.436517000000,0.979460000000,0.008872000000
+545,0.451584000000,0.982200000000,0.007918000000
+546,0.466782000000,0.984520000000,0.007030000000
+547,0.482147000000,0.986520000000,0.006223000000
+548,0.497738000000,0.988320000000,0.005453000000
+549,0.513606000000,0.990020000000,0.004714000000
+550,0.529826000000,0.991761000000,0.003988000000
+551,0.546440000000,0.993530000000,0.003289000000
+552,0.563426000000,0.995230000000,0.002646000000
+553,0.580726000000,0.996770000000,0.002063000000
+554,0.598290000000,0.998090000000,0.001533000000
+555,0.616053000000,0.999110000000,0.001091000000
+556,0.633948000000,0.999770000000,0.000711000000
+557,0.651901000000,1.000000000000,0.000407000000
+558,0.669824000000,0.999710000000,0.000184000000
+559,0.687632000000,0.998850000000,0.000047000000
+560,0.705224000000,0.997340000000,0.000000000000
+561,0.722773000000,0.995260000000,0.000000000000
+562,0.740483000000,0.992740000000,0.000000000000
+563,0.758273000000,0.989750000000,0.000000000000
+564,0.776083000000,0.986300000000,0.000000000000
+565,0.793832000000,0.982380000000,0.000000000000
+566,0.811436000000,0.977980000000,0.000000000000
+567,0.828822000000,0.973110000000,0.000000000000
+568,0.845879000000,0.967740000000,0.000000000000
+569,0.862525000000,0.961890000000,0.000000000000
+570,0.878655000000,0.955552000000,0.000000000000
+571,0.894208000000,0.948601000000,0.000000000000
+572,0.909206000000,0.940981000000,0.000000000000
+573,0.923672000000,0.932798000000,0.000000000000
+574,0.937638000000,0.924158000000,0.000000000000
+575,0.951162000000,0.915175000000,0.000000000000
+576,0.964283000000,0.905954000000,0.000000000000
+577,0.977068000000,0.896608000000,0.000000000000
+578,0.989590000000,0.887249000000,0.000000000000
+579,1.001910000000,0.877986000000,0.000000000000
+580,1.014160000000,0.868934000000,0.000000000000
+581,1.026500000000,0.860164000000,0.000000000000
+582,1.038800000000,0.851519000000,0.000000000000
+583,1.051000000000,0.842963000000,0.000000000000
+584,1.062900000000,0.834393000000,0.000000000000
+585,1.074300000000,0.825623000000,0.000000000000
+586,1.085200000000,0.816764000000,0.000000000000
+587,1.095200000000,0.807544000000,0.000000000000
+588,1.104200000000,0.797947000000,0.000000000000
+589,1.112000000000,0.787893000000,0.000000000000
+590,1.118520000000,0.777405000000,0.000000000000
+591,1.123800000000,0.766490000000,0.000000000000
+592,1.128000000000,0.755309000000,0.000000000000
+593,1.131100000000,0.743845000000,0.000000000000
+594,1.133200000000,0.732190000000,0.000000000000
+595,1.134300000000,0.720353000000,0.000000000000
+596,1.134300000000,0.708281000000,0.000000000000
+597,1.133300000000,0.696055000000,0.000000000000
+598,1.131200000000,0.683621000000,0.000000000000
+599,1.128100000000,0.671048000000,0.000000000000
+600,1.123990000000,0.658341000000,0.000000000000
+601,1.118900000000,0.645545000000,0.000000000000
+602,1.112900000000,0.632718000000,0.000000000000
+603,1.105900000000,0.619815000000,0.000000000000
+604,1.098000000000,0.606887000000,0.000000000000
+605,1.089100000000,0.593878000000,0.000000000000
+606,1.079200000000,0.580781000000,0.000000000000
+607,1.068400000000,0.567653000000,0.000000000000
+608,1.056700000000,0.554490000000,0.000000000000
+609,1.044000000000,0.541228000000,0.000000000000
+610,1.030480000000,0.527963000000,0.000000000000
+611,1.016000000000,0.514634000000,0.000000000000
+612,1.000800000000,0.501363000000,0.000000000000
+613,0.984790000000,0.488124000000,0.000000000000
+614,0.968080000000,0.474935000000,0.000000000000
+615,0.950740000000,0.461834000000,0.000000000000
+616,0.932800000000,0.448823000000,0.000000000000
+617,0.914340000000,0.435917000000,0.000000000000
+618,0.895390000000,0.423153000000,0.000000000000
+619,0.876030000000,0.410526000000,0.000000000000
+620,0.856297000000,0.398057000000,0.000000000000
+621,0.836350000000,0.385835000000,0.000000000000
+622,0.816290000000,0.373951000000,0.000000000000
+623,0.796050000000,0.362311000000,0.000000000000
+624,0.775610000000,0.350863000000,0.000000000000
+625,0.754930000000,0.339554000000,0.000000000000
+626,0.733990000000,0.328309000000,0.000000000000
+627,0.712780000000,0.317118000000,0.000000000000
+628,0.691290000000,0.305936000000,0.000000000000
+629,0.669520000000,0.294737000000,0.000000000000
+630,0.647467000000,0.283493000000,0.000000000000
+631,0.625110000000,0.272222000000,0.000000000000
+632,0.602520000000,0.260990000000,0.000000000000
+633,0.579890000000,0.249877000000,0.000000000000
+634,0.557370000000,0.238946000000,0.000000000000
+635,0.535110000000,0.228254000000,0.000000000000
+636,0.513240000000,0.217853000000,0.000000000000
+637,0.491860000000,0.207780000000,0.000000000000
+638,0.471080000000,0.198072000000,0.000000000000
+639,0.450960000000,0.188748000000,0.000000000000
+640,0.431567000000,0.179828000000,0.000000000000
+641,0.412870000000,0.171285000000,0.000000000000
+642,0.394750000000,0.163059000000,0.000000000000
+643,0.377210000000,0.155151000000,0.000000000000
+644,0.360190000000,0.147535000000,0.000000000000
+645,0.343690000000,0.140211000000,0.000000000000
+646,0.327690000000,0.133170000000,0.000000000000
+647,0.312170000000,0.126400000000,0.000000000000
+648,0.297110000000,0.119892000000,0.000000000000
+649,0.282500000000,0.113640000000,0.000000000000
+650,0.268329000000,0.107633000000,0.000000000000
+651,0.254590000000,0.101870000000,0.000000000000
+652,0.241300000000,0.096347000000,0.000000000000
+653,0.228480000000,0.091063000000,0.000000000000
+654,0.216140000000,0.086010000000,0.000000000000
+655,0.204300000000,0.081187000000,0.000000000000
+656,0.192950000000,0.076583000000,0.000000000000
+657,0.182110000000,0.072198000000,0.000000000000
+658,0.171770000000,0.068024000000,0.000000000000
+659,0.161920000000,0.064052000000,0.000000000000
+660,0.152568000000,0.060281000000,0.000000000000
+661,0.143670000000,0.056697000000,0.000000000000
+662,0.135200000000,0.053292000000,0.000000000000
+663,0.127130000000,0.050059000000,0.000000000000
+664,0.119480000000,0.046998000000,0.000000000000
+665,0.112210000000,0.044096000000,0.000000000000
+666,0.105310000000,0.041345000000,0.000000000000
+667,0.098786000000,0.038750700000,0.000000000000
+668,0.092610000000,0.036297800000,0.000000000000
+669,0.086773000000,0.033983200000,0.000000000000
+670,0.081260600000,0.031800400000,0.000000000000
+671,0.076048000000,0.029739500000,0.000000000000
+672,0.071114000000,0.027791800000,0.000000000000
+673,0.066454000000,0.025955100000,0.000000000000
+674,0.062062000000,0.024226300000,0.000000000000
+675,0.057930000000,0.022601700000,0.000000000000
+676,0.054050000000,0.021077900000,0.000000000000
+677,0.050412000000,0.019650500000,0.000000000000
+678,0.047006000000,0.018315300000,0.000000000000
+679,0.043823000000,0.017068600000,0.000000000000
+680,0.040850800000,0.015905100000,0.000000000000
+681,0.038072000000,0.014818300000,0.000000000000
+682,0.035468000000,0.013800800000,0.000000000000
+683,0.033031000000,0.012849500000,0.000000000000
+684,0.030753000000,0.011960700000,0.000000000000
+685,0.028623000000,0.011130300000,0.000000000000
+686,0.026635000000,0.010355500000,0.000000000000
+687,0.024781000000,0.009633200000,0.000000000000
+688,0.023052000000,0.008959900000,0.000000000000
+689,0.021441000000,0.008332400000,0.000000000000
+690,0.019941300000,0.007748800000,0.000000000000
+691,0.018544000000,0.007204600000,0.000000000000
+692,0.017241000000,0.006697500000,0.000000000000
+693,0.016027000000,0.006225100000,0.000000000000
+694,0.014896000000,0.005785000000,0.000000000000
+695,0.013842000000,0.005375100000,0.000000000000
+696,0.012862000000,0.004994100000,0.000000000000
+697,0.011949000000,0.004639200000,0.000000000000
+698,0.011100000000,0.004309300000,0.000000000000
+699,0.010311000000,0.004002800000,0.000000000000
+700,0.009576880000,0.003717740000,0.000000000000
+701,0.008894000000,0.003452620000,0.000000000000
+702,0.008258100000,0.003205830000,0.000000000000
+703,0.007666400000,0.002976230000,0.000000000000
+704,0.007116300000,0.002762810000,0.000000000000
+705,0.006605200000,0.002564560000,0.000000000000
+706,0.006130600000,0.002380480000,0.000000000000
+707,0.005690300000,0.002209710000,0.000000000000
+708,0.005281900000,0.002051320000,0.000000000000
+709,0.004903300000,0.001904490000,0.000000000000
+710,0.004552630000,0.001768470000,0.000000000000
+711,0.004227500000,0.001642360000,0.000000000000
+712,0.003925800000,0.001525350000,0.000000000000
+713,0.003645700000,0.001416720000,0.000000000000
+714,0.003385900000,0.001315950000,0.000000000000
+715,0.003144700000,0.001222390000,0.000000000000
+716,0.002920800000,0.001135550000,0.000000000000
+717,0.002713000000,0.001054940000,0.000000000000
+718,0.002520200000,0.000980140000,0.000000000000
+719,0.002341100000,0.000910660000,0.000000000000
+720,0.002174960000,0.000846190000,0.000000000000
+721,0.002020600000,0.000786290000,0.000000000000
+722,0.001877300000,0.000730680000,0.000000000000
+723,0.001744100000,0.000678990000,0.000000000000
+724,0.001620500000,0.000631010000,0.000000000000
+725,0.001505700000,0.000586440000,0.000000000000
+726,0.001399200000,0.000545110000,0.000000000000
+727,0.001300400000,0.000506720000,0.000000000000
+728,0.001208700000,0.000471110000,0.000000000000
+729,0.001123600000,0.000438050000,0.000000000000
+730,0.001044760000,0.000407410000,0.000000000000
+731,0.000971560000,0.000378962000,0.000000000000
+732,0.000903600000,0.000352543000,0.000000000000
+733,0.000840480000,0.000328001000,0.000000000000
+734,0.000781870000,0.000305208000,0.000000000000
+735,0.000727450000,0.000284041000,0.000000000000
+736,0.000676900000,0.000264375000,0.000000000000
+737,0.000629960000,0.000246109000,0.000000000000
+738,0.000586370000,0.000229143000,0.000000000000
+739,0.000545870000,0.000213376000,0.000000000000
+740,0.000508258000,0.000198730000,0.000000000000
+741,0.000473300000,0.000185115000,0.000000000000
+742,0.000440800000,0.000172454000,0.000000000000
+743,0.000410580000,0.000160678000,0.000000000000
+744,0.000382490000,0.000149730000,0.000000000000
+745,0.000356380000,0.000139550000,0.000000000000
+746,0.000332110000,0.000130086000,0.000000000000
+747,0.000309550000,0.000121290000,0.000000000000
+748,0.000288580000,0.000113106000,0.000000000000
+749,0.000269090000,0.000105501000,0.000000000000
+750,0.000250969000,0.000098428000,0.000000000000
+751,0.000234130000,0.000091853000,0.000000000000
+752,0.000218470000,0.000085738000,0.000000000000
+753,0.000203910000,0.000080048000,0.000000000000
+754,0.000190350000,0.000074751000,0.000000000000
+755,0.000177730000,0.000069819000,0.000000000000
+756,0.000165970000,0.000065222000,0.000000000000
+757,0.000155020000,0.000060939000,0.000000000000
+758,0.000144800000,0.000056942000,0.000000000000
+759,0.000135280000,0.000053217000,0.000000000000
+760,0.000126390000,0.000049737000,0.000000000000
+761,0.000118100000,0.000046491000,0.000000000000
+762,0.000110370000,0.000043464000,0.000000000000
+763,0.000103150000,0.000040635000,0.000000000000
+764,0.000096427000,0.000038000000,0.000000000000
+765,0.000090151000,0.000035540500,0.000000000000
+766,0.000084294000,0.000033244800,0.000000000000
+767,0.000078830000,0.000031100600,0.000000000000
+768,0.000073729000,0.000029099000,0.000000000000
+769,0.000068969000,0.000027230700,0.000000000000
+770,0.000064525800,0.000025486000,0.000000000000
+771,0.000060376000,0.000023856100,0.000000000000
+772,0.000056500000,0.000022333200,0.000000000000
+773,0.000052880000,0.000020910400,0.000000000000
+774,0.000049498000,0.000019580800,0.000000000000
+775,0.000046339000,0.000018338400,0.000000000000
+776,0.000043389000,0.000017177700,0.000000000000
+777,0.000040634000,0.000016093400,0.000000000000
+778,0.000038060000,0.000015080000,0.000000000000
+779,0.000035657000,0.000014133600,0.000000000000
+780,0.000033411700,0.000013249000,0.000000000000
+781,0.000031315000,0.000012422600,0.000000000000
+782,0.000029355000,0.000011649900,0.000000000000
+783,0.000027524000,0.000010927700,0.000000000000
+784,0.000025811000,0.000010251900,0.000000000000
+785,0.000024209000,0.000009619600,0.000000000000
+786,0.000022711000,0.000009028100,0.000000000000
+787,0.000021308000,0.000008474000,0.000000000000
+788,0.000019994000,0.000007954800,0.000000000000
+789,0.000018764000,0.000007468600,0.000000000000
+790,0.000017611500,0.000007012800,0.000000000000
+791,0.000016532000,0.000006585800,0.000000000000
+792,0.000015521000,0.000006185700,0.000000000000
+793,0.000014574000,0.000005810700,0.000000000000
+794,0.000013686000,0.000005459000,0.000000000000
+795,0.000012855000,0.000005129800,0.000000000000
+796,0.000012075000,0.000004820600,0.000000000000
+797,0.000011345000,0.000004531200,0.000000000000
+798,0.000010659000,0.000004259100,0.000000000000
+799,0.000010017000,0.000004004200,0.000000000000
+800,0.000009413630,0.000003764730,0.000000000000
+801,0.000008847900,0.000003539950,0.000000000000
+802,0.000008317100,0.000003329140,0.000000000000
+803,0.000007819000,0.000003131150,0.000000000000
+804,0.000007351600,0.000002945290,0.000000000000
+805,0.000006913000,0.000002770810,0.000000000000
+806,0.000006501500,0.000002607050,0.000000000000
+807,0.000006115300,0.000002453290,0.000000000000
+808,0.000005752900,0.000002308940,0.000000000000
+809,0.000005412700,0.000002173380,0.000000000000
+810,0.000005093470,0.000002046130,0.000000000000
+811,0.000004793800,0.000001926620,0.000000000000
+812,0.000004512500,0.000001814400,0.000000000000
+813,0.000004248300,0.000001708950,0.000000000000
+814,0.000004000200,0.000001609880,0.000000000000
+815,0.000003767100,0.000001516770,0.000000000000
+816,0.000003548000,0.000001429210,0.000000000000
+817,0.000003342100,0.000001346860,0.000000000000
+818,0.000003148500,0.000001269450,0.000000000000
+819,0.000002966500,0.000001196620,0.000000000000
+820,0.000002795310,0.000001128090,0.000000000000
+821,0.000002634500,0.000001063680,0.000000000000
+822,0.000002483400,0.000001003130,0.000000000000
+823,0.000002341400,0.000000946220,0.000000000000
+824,0.000002207800,0.000000892630,0.000000000000
+825,0.000002082000,0.000000842160,0.000000000000
+826,0.000001963600,0.000000794640,0.000000000000
+827,0.000001851900,0.000000749780,0.000000000000
+828,0.000001746500,0.000000707440,0.000000000000
+829,0.000001647100,0.000000667480,0.000000000000
+830,0.000001553140,0.000000629700,0.000000000000
\ No newline at end of file diff --git a/geometry.py b/geometry.py index 8af8d49..8dee8c9 100644 --- a/geometry.py +++ b/geometry.py @@ -106,25 +106,17 @@ class Geometry(object): else: self.material2_index[i] = -1 - surface1 = np.concatenate([solid.surface1 for solid in self.solids]) - surface2 = np.concatenate([solid.surface2 for solid in self.solids]) + surface = np.concatenate([solid.surface for solid in self.solids]) - self.surfaces = list(np.unique(np.concatenate([surface1, surface2]))) + self.surfaces = list(np.unique(surface)) - self.surface1_index = np.empty(len(self.mesh), dtype=np.int32) - self.surface2_index = np.empty(len(self.mesh), dtype=np.int32) + self.surface_index = np.empty(len(self.mesh), dtype=np.int32) - for i, surface in enumerate(surface1[reorder]): + for i, surface in enumerate(surface[reorder]): if surface is not None: - self.surface1_index[i] = self.surfaces.index(surface) + self.surface_index[i] = self.surfaces.index(surface) else: - self.surface1_index[i] = -1 - - for i, surface in enumerate(surface2[reorder]): - if surface is not None: - self.surface2_index[i] = self.surfaces.index(surface) - else: - self.surface2_index[i] = -1 + self.surface_index[i] = -1 self.colors = np.concatenate([solid.color for solid in self.solids]) @@ -210,6 +202,10 @@ class Geometry(object): of the texture references. """ + set_wavelength_range = module.get_function('set_wavelength_range') + + set_wavelength_range(wavelengths[0], wavelengths[-1], wavelengths[1]-wavelengths[0], block=(1,1,1), grid=(1,1)) + set_material = module.get_function('set_material') for i, material in enumerate(self.materials): if material is None: @@ -223,7 +219,7 @@ class Geometry(object): material.absorption_length_gpu = cuda.to_device(absorption_length) material.scattering_length_gpu = cuda.to_device(scattering_length) - set_material(i, material.refractive_index_gpu, material.absorption_length_gpu, material.scattering_length_gpu, block=(1,1,1), grid=(1,1)) + set_material(np.int32(i), material.refractive_index_gpu, material.absorption_length_gpu, material.scattering_length_gpu, block=(1,1,1), grid=(1,1)) set_surface = module.get_function('set_surface') for i, surface in enumerate(self.surfaces): @@ -232,33 +228,29 @@ class Geometry(object): absorption = np.interp(wavelengths, surface.absorption[:,0], surface.absorption[:,1]).astype(np.float32) reflection_diffuse = np.interp(wavelengths, surface.reflection_diffuse[:,0], surface.reflection_diffuse[:,1]).astype(np.float32) - reflection_specular = np.interp(wavelengths, surface.refelection_specular[:,0], surface.reflection_specular[:,1]).astype(np.float32) + reflection_specular = np.interp(wavelengths, surface.reflection_specular[:,0], surface.reflection_specular[:,1]).astype(np.float32) surface.absorption_gpu = cuda.to_device(absorption) surface.reflection_diffuse_gpu = cuda.to_device(reflection_diffuse) surface.reflection_specular_gpu = cuda.to_device(reflection_specular) - set_surface(i, surface.absorption_gpu, surface.reflection_diffuse_gpu, surface.reflection_specular_gpu) + set_surface(np.int32(i), surface.absorption_gpu, surface.reflection_diffuse_gpu, surface.reflection_specular_gpu, block=(1,1,1), grid=(1,1)) - material1_index_tex = module.get_texref('material1_index') - material2_index_tex = module.get_texref('material2_index') - surface1_index_tex = module.get_texref('surface1_index') - surface2_index_tex = module.get_texref('surface2_index') + material1_index_tex = module.get_texref('material1_lookup') + material2_index_tex = module.get_texref('material2_lookup') + surface_index_tex = module.get_texref('surface_lookup') material1_index_gpu = cuda.to_device(self.material1_index) material2_index_gpu = cuda.to_device(self.material2_index) - surface1_index_gpu = cuda.to_device(self.surface1_index) - surface2_index_gpu = cuda.to_device(self.surface2_index) + surface_index_gpu = cuda.to_device(self.surface_index) material1_index_tex.set_address(material1_index_gpu, self.material1_index.nbytes) material2_index_tex.set_address(material2_index_gpu, self.material2_index.nbytes) - surface1_index_tex.set_address(surface1_index_gpu, self.surface1_index.nbytes) - surface2_index_tex.set_address(surface2_index_gpu, self.surface2_index.nbytes) + surface_index_tex.set_address(surface_index_gpu, self.surface_index.nbytes) material1_index_tex.set_format(cuda.array_format.SIGNED_INT32, 1) material2_index_tex.set_format(cuda.array_format.SIGNED_INT32, 1) - surface1_index_tex.set_format(cuda.array_format.SIGNED_INT32, 1) - surface2_index_tex.set_format(cuda.array_format.SIGNED_INT32, 1) + surface_index_tex.set_format(cuda.array_format.SIGNED_INT32, 1) vertices = np.empty(len(self.mesh.vertices), dtype=gpuarray.vec.float4) vertices['x'] = self.mesh.vertices[:,0] @@ -274,7 +266,7 @@ class Geometry(object): if color: triangles['w'] = self.colors else: - triangles['w'] = (self.material1_index << 24) | (self.material2_index << 16) | (self.surface1_index << 8) | self.surface2_index + triangles['w'] = (self.material1_index << 24) | (self.material2_index << 16) | (self.surface_index << 8) lower_bounds = np.empty(self.lower_bounds.shape[0], dtype=gpuarray.vec.float4) lower_bounds['x'] = self.lower_bounds[:,0] diff --git a/materials.py b/materials.py index dc28a8b..273a4bb 100644 --- a/materials.py +++ b/materials.py @@ -1,6 +1,6 @@ import numpy as np -wavelengths = np.arange(200, 810, 20) +wavelengths = np.arange(200, 810, 20).astype(np.float32) class Material(object): """Material optical properties.""" @@ -68,10 +68,10 @@ r1408_photocathode = Surface('photocathode') r1408_photocathode.absorption = \ np.array([(60.0, 0.0), (200.0, 0.0), (230.0, 0.0), (231.0, 0.0), (232.0, 0.0), (233.0, 0.0), (234.0, 0.0), (235.0, 0.0), (236.0, 0.0), (237.0, 0.0), (238.0, 0.0), (239.0, 0.0), (240.0, 0.0), (241.0, 0.00049089999999999995), (242.0, 0.0009817999999999999), (243.0, 0.0014729999999999999), (244.0, 0.001964), (245.0, 0.002454), (246.0, 0.0029450000000000001), (247.0, 0.0034359999999999998), (248.0, 0.0039269999999999999), (249.0, 0.0044180000000000001), (250.0, 0.0049090000000000002), (251.0, 0.0055100000000000001), (252.0, 0.0061110000000000001), (253.0, 0.0067130000000000002), (254.0, 0.0073140000000000002), (255.0, 0.0079150000000000002), (256.0, 0.0085159999999999993), (257.0, 0.0091170000000000001), (258.0, 0.0097190000000000002), (259.0, 0.010319999999999999), (260.0, 0.010919999999999999), (261.0, 0.01269), (262.0, 0.014460000000000001), (263.0, 0.016219999999999998), (264.0, 0.017989999999999999), (265.0, 0.01976), (266.0, 0.021530000000000001), (267.0, 0.023290000000000002), (268.0, 0.025059999999999999), (269.0, 0.02683), (270.0, 0.0286), (271.0, 0.031690000000000003), (272.0, 0.034779999999999998), (273.0, 0.037879999999999997), (274.0, 0.040969999999999999), (275.0, 0.044060000000000002), (276.0, 0.047160000000000001), (277.0, 0.050250000000000003), (278.0, 0.053339999999999999), (279.0, 0.056439999999999997), (280.0, 0.05953), (281.0, 0.062950000000000006), (282.0, 0.066360000000000002), (283.0, 0.069779999999999995), (284.0, 0.073200000000000001), (285.0, 0.076609999999999998), (286.0, 0.080030000000000004), (287.0, 0.08344), (288.0, 0.086860000000000007), (289.0, 0.090279999999999999), (290.0, 0.093689999999999996), (291.0, 0.096799999999999997), (292.0, 0.099909999999999999), (293.0, 0.10299999999999999), (294.0, 0.1061), (295.0, 0.10920000000000001), (296.0, 0.1123), (297.0, 0.1154), (298.0, 0.11849999999999999), (299.0, 0.1217), (300.0, 0.12479999999999999), (301.0, 0.1275), (302.0, 0.13020000000000001), (303.0, 0.13300000000000001), (304.0, 0.13569999999999999), (305.0, 0.13850000000000001), (306.0, 0.14119999999999999), (307.0, 0.1439), (308.0, 0.1467), (309.0, 0.14940000000000001), (310.0, 0.1522), (311.0, 0.15409999999999999), (312.0, 0.15609999999999999), (313.0, 0.158), (314.0, 0.16), (315.0, 0.16189999999999999), (316.0, 0.16389999999999999), (317.0, 0.1658), (318.0, 0.1678), (319.0, 0.16969999999999999), (320.0, 0.17169999999999999), (321.0, 0.17249999999999999), (322.0, 0.1734), (323.0, 0.17419999999999999), (324.0, 0.17510000000000001), (325.0, 0.1759), (326.0, 0.1767), (327.0, 0.17760000000000001), (328.0, 0.1784), (329.0, 0.17929999999999999), (330.0, 0.18010000000000001), (331.0, 0.18099999999999999), (332.0, 0.18190000000000001), (333.0, 0.1827), (334.0, 0.18360000000000001), (335.0, 0.18440000000000001), (336.0, 0.18529999999999999), (337.0, 0.1862), (338.0, 0.187), (339.0, 0.18790000000000001), (340.0, 0.1888), (341.0, 0.189), (342.0, 0.1893), (343.0, 0.18959999999999999), (344.0, 0.1898), (345.0, 0.19009999999999999), (346.0, 0.1903), (347.0, 0.19059999999999999), (348.0, 0.19089999999999999), (349.0, 0.19109999999999999), (350.0, 0.19139999999999999), (351.0, 0.19170000000000001), (352.0, 0.19209999999999999), (353.0, 0.19239999999999999), (354.0, 0.1928), (355.0, 0.19309999999999999), (356.0, 0.19350000000000001), (357.0, 0.1938), (358.0, 0.19420000000000001), (359.0, 0.19450000000000001), (360.0, 0.19489999999999999), (361.0, 0.19520000000000001), (362.0, 0.19539999999999999), (363.0, 0.19570000000000001), (364.0, 0.19600000000000001), (365.0, 0.1963), (366.0, 0.1966), (367.0, 0.19689999999999999), (368.0, 0.19719999999999999), (369.0, 0.19750000000000001), (370.0, 0.1978), (371.0, 0.1978), (372.0, 0.1978), (373.0, 0.1978), (374.0, 0.19789999999999999), (375.0, 0.19789999999999999), (376.0, 0.19789999999999999), (377.0, 0.19789999999999999), (378.0, 0.19800000000000001), (379.0, 0.19800000000000001), (380.0, 0.19800000000000001), (381.0, 0.1983), (382.0, 0.19850000000000001), (383.0, 0.1988), (384.0, 0.19900000000000001), (385.0, 0.1993), (386.0, 0.1996), (387.0, 0.19980000000000001), (388.0, 0.2001), (389.0, 0.20030000000000001), (390.0, 0.2006), (391.0, 0.20080000000000001), (392.0, 0.20100000000000001), (393.0, 0.20130000000000001), (394.0, 0.20150000000000001), (395.0, 0.20169999999999999), (396.0, 0.2019), (397.0, 0.20219999999999999), (398.0, 0.2024), (399.0, 0.2026), (400.0, 0.2029), (401.0, 0.20330000000000001), (402.0, 0.20380000000000001), (403.0, 0.20419999999999999), (404.0, 0.20469999999999999), (405.0, 0.20519999999999999), (406.0, 0.2056), (407.0, 0.20610000000000001), (408.0, 0.20649999999999999), (409.0, 0.20699999999999999), (410.0, 0.20749999999999999), (411.0, 0.20780000000000001), (412.0, 0.2082), (413.0, 0.20849999999999999), (414.0, 0.2089), (415.0, 0.2092), (416.0, 0.20960000000000001), (417.0, 0.2099), (418.0, 0.21029999999999999), (419.0, 0.21060000000000001), (420.0, 0.21099999999999999), (421.0, 0.21129999999999999), (422.0, 0.21160000000000001), (423.0, 0.21190000000000001), (424.0, 0.21229999999999999), (425.0, 0.21260000000000001), (426.0, 0.21290000000000001), (427.0, 0.2132), (428.0, 0.2135), (429.0, 0.21390000000000001), (430.0, 0.2142), (431.0, 0.2142), (432.0, 0.21429999999999999), (433.0, 0.21440000000000001), (434.0, 0.21440000000000001), (435.0, 0.2145), (436.0, 0.21460000000000001), (437.0, 0.21460000000000001), (438.0, 0.2147), (439.0, 0.21479999999999999), (440.0, 0.21479999999999999), (441.0, 0.2147), (442.0, 0.2145), (443.0, 0.21440000000000001), (444.0, 0.2142), (445.0, 0.21410000000000001), (446.0, 0.21390000000000001), (447.0, 0.21379999999999999), (448.0, 0.21360000000000001), (449.0, 0.2135), (450.0, 0.21329999999999999), (451.0, 0.21310000000000001), (452.0, 0.21279999999999999), (453.0, 0.21260000000000001), (454.0, 0.21240000000000001), (455.0, 0.21210000000000001), (456.0, 0.21190000000000001), (457.0, 0.2117), (458.0, 0.2114), (459.0, 0.2112), (460.0, 0.21099999999999999), (461.0, 0.21029999999999999), (462.0, 0.2097), (463.0, 0.20899999999999999), (464.0, 0.2084), (465.0, 0.2077), (466.0, 0.20710000000000001), (467.0, 0.2064), (468.0, 0.20580000000000001), (469.0, 0.2051), (470.0, 0.20449999999999999), (471.0, 0.20380000000000001), (472.0, 0.20319999999999999), (473.0, 0.20250000000000001), (474.0, 0.20180000000000001), (475.0, 0.20119999999999999), (476.0, 0.20050000000000001), (477.0, 0.19989999999999999), (478.0, 0.19919999999999999), (479.0, 0.1986), (480.0, 0.19789999999999999), (481.0, 0.19750000000000001), (482.0, 0.1971), (483.0, 0.19670000000000001), (484.0, 0.19620000000000001), (485.0, 0.1958), (486.0, 0.19539999999999999), (487.0, 0.19500000000000001), (488.0, 0.1946), (489.0, 0.19409999999999999), (490.0, 0.19370000000000001), (491.0, 0.1933), (492.0, 0.19289999999999999), (493.0, 0.1925), (494.0, 0.19209999999999999), (495.0, 0.19170000000000001), (496.0, 0.1913), (497.0, 0.1908), (498.0, 0.19040000000000001), (499.0, 0.19), (500.0, 0.18959999999999999), (501.0, 0.18920000000000001), (502.0, 0.18870000000000001), (503.0, 0.1883), (504.0, 0.18790000000000001), (505.0, 0.18740000000000001), (506.0, 0.187), (507.0, 0.18659999999999999), (508.0, 0.18609999999999999), (509.0, 0.1857), (510.0, 0.1852), (511.0, 0.1835), (512.0, 0.18179999999999999), (513.0, 0.18010000000000001), (514.0, 0.1784), (515.0, 0.17660000000000001), (516.0, 0.1749), (517.0, 0.17319999999999999), (518.0, 0.17150000000000001), (519.0, 0.16980000000000001), (520.0, 0.1681), (521.0, 0.1651), (522.0, 0.16220000000000001), (523.0, 0.1593), (524.0, 0.15640000000000001), (525.0, 0.1535), (526.0, 0.15060000000000001), (527.0, 0.1477), (528.0, 0.14480000000000001), (529.0, 0.1419), (530.0, 0.13900000000000001), (531.0, 0.1366), (532.0, 0.13420000000000001), (533.0, 0.1318), (534.0, 0.12939999999999999), (535.0, 0.127), (536.0, 0.1246), (537.0, 0.1222), (538.0, 0.1198), (539.0, 0.1174), (540.0, 0.115), (541.0, 0.1143), (542.0, 0.11360000000000001), (543.0, 0.1129), (544.0, 0.11219999999999999), (545.0, 0.1115), (546.0, 0.11070000000000001), (547.0, 0.11), (548.0, 0.10929999999999999), (549.0, 0.1086), (550.0, 0.1079), (551.0, 0.1062), (552.0, 0.1046), (553.0, 0.10290000000000001), (554.0, 0.1013), (555.0, 0.099650000000000002), (556.0, 0.09801), (557.0, 0.096360000000000001), (558.0, 0.094719999999999999), (559.0, 0.093079999999999996), (560.0, 0.091429999999999997), (561.0, 0.0906), (562.0, 0.089770000000000003), (563.0, 0.088950000000000001), (564.0, 0.088120000000000004), (565.0, 0.087290000000000006), (566.0, 0.086459999999999995), (567.0, 0.085629999999999998), (568.0, 0.0848), (569.0, 0.083979999999999999), (570.0, 0.083150000000000002), (571.0, 0.082239999999999994), (572.0, 0.08133), (573.0, 0.080409999999999995), (574.0, 0.079500000000000001), (575.0, 0.078589999999999993), (576.0, 0.077679999999999999), (577.0, 0.076770000000000005), (578.0, 0.075859999999999997), (579.0, 0.074950000000000003), (580.0, 0.074039999999999995), (581.0, 0.072969999999999993), (582.0, 0.071900000000000006), (583.0, 0.070830000000000004), (584.0, 0.069760000000000003), (585.0, 0.068690000000000001), (586.0, 0.06762), (587.0, 0.066549999999999998), (588.0, 0.065479999999999997), (589.0, 0.064409999999999995), (590.0, 0.063339999999999994), (591.0, 0.062230000000000001), (592.0, 0.061109999999999998), (593.0, 0.059990000000000002), (594.0, 0.058869999999999999), (595.0, 0.057750000000000003), (596.0, 0.05663), (597.0, 0.05552), (598.0, 0.054399999999999997), (599.0, 0.053280000000000001), (600.0, 0.052159999999999998), (601.0, 0.05092), (602.0, 0.049669999999999999), (603.0, 0.048430000000000001), (604.0, 0.047190000000000003), (605.0, 0.045940000000000002), (606.0, 0.044699999999999997), (607.0, 0.043450000000000003), (608.0, 0.042209999999999998), (609.0, 0.040969999999999999), (610.0, 0.039719999999999998), (611.0, 0.03884), (612.0, 0.037949999999999998), (613.0, 0.037069999999999999), (614.0, 0.036179999999999997), (615.0, 0.035299999999999998), (616.0, 0.034410000000000003), (617.0, 0.033529999999999997), (618.0, 0.032640000000000002), (619.0, 0.031759999999999997), (620.0, 0.030870000000000002), (621.0, 0.029960000000000001), (622.0, 0.029059999999999999), (623.0, 0.028150000000000001), (624.0, 0.02724), (625.0, 0.026339999999999999), (626.0, 0.025430000000000001), (627.0, 0.02452), (628.0, 0.023609999999999999), (629.0, 0.022710000000000001), (630.0, 0.0218), (631.0, 0.021090000000000001), (632.0, 0.020379999999999999), (633.0, 0.01967), (634.0, 0.018970000000000001), (635.0, 0.018259999999999998), (636.0, 0.01755), (637.0, 0.016840000000000001), (638.0, 0.016129999999999999), (639.0, 0.01542), (640.0, 0.014710000000000001), (641.0, 0.014120000000000001), (642.0, 0.01353), (643.0, 0.012930000000000001), (644.0, 0.01234), (645.0, 0.01175), (646.0, 0.01115), (647.0, 0.01056), (648.0, 0.0099670000000000002), (649.0, 0.0093740000000000004), (650.0, 0.0087810000000000006), (651.0, 0.0084159999999999999), (652.0, 0.0080520000000000001), (653.0, 0.0076870000000000003), (654.0, 0.0073229999999999996), (655.0, 0.0069579999999999998), (656.0, 0.006594), (657.0, 0.0062290000000000002), (658.0, 0.0058650000000000004), (659.0, 0.0054999999999999997), (660.0, 0.0051359999999999999), (661.0, 0.0048989999999999997), (662.0, 0.0046629999999999996), (663.0, 0.0044270000000000004), (664.0, 0.0041910000000000003), (665.0, 0.0039550000000000002), (666.0, 0.003718), (667.0, 0.0034819999999999999), (668.0, 0.0032460000000000002), (669.0, 0.0030100000000000001), (670.0, 0.002774), (671.0, 0.00265), (672.0, 0.0025270000000000002), (673.0, 0.0024039999999999999), (674.0, 0.002281), (675.0, 0.0021570000000000001), (676.0, 0.0020339999999999998), (677.0, 0.0019109999999999999), (678.0, 0.0017880000000000001), (679.0, 0.001665), (680.0, 0.0015410000000000001), (681.0, 0.001459), (682.0, 0.001377), (683.0, 0.0012949999999999999), (684.0, 0.0012130000000000001), (685.0, 0.0011310000000000001), (686.0, 0.001049), (687.0, 0.00096650000000000002), (688.0, 0.00088440000000000003), (689.0, 0.00080230000000000004), (690.0, 0.00072009999999999999), (691.0, 0.00069439999999999997), (692.0, 0.00066859999999999999), (693.0, 0.00064289999999999996), (694.0, 0.00061709999999999998), (695.0, 0.00059139999999999996), (696.0, 0.00056559999999999998), (697.0, 0.00053989999999999995), (698.0, 0.00051409999999999997), (699.0, 0.00048840000000000005), (700.0, 0.00046260000000000002), (800.0, 0.0)]) -r1408_photocathode.specular_reflection = \ +r1408_photocathode.reflection_specular = \ np.array(zip(wavelengths, np.zeros(len(wavelengths)))) -r1408_photocathode.diffuse_reflection = \ +r1408_photocathode.reflection_diffuse = \ np.array(zip(wavelengths, np.zeros(len(wavelengths)))) aluminum = Surface('aluminum') @@ -79,8 +79,8 @@ aluminum = Surface('aluminum') aluminum.absorption = \ np.array(zip(wavelengths, np.zeros(len(wavelengths)))) -aluminum.specular_reflection = \ +aluminum.reflection_specular = \ np.array(zip(wavelengths, np.ones(len(wavelengths)))) -aluminum.diffuse_reflection = \ +aluminum.reflection_diffuse = \ np.array(zip(wavelengths, np.zeros(len(wavelengths)))) diff --git a/models/Colbert_HighRes_Brow.STL b/models/Colbert_HighRes_Brow.STL Binary files differnew file mode 100644 index 0000000..784c877 --- /dev/null +++ b/models/Colbert_HighRes_Brow.STL diff --git a/models/Colbert_HighRes_Smile.STL b/models/Colbert_HighRes_Smile.STL Binary files differnew file mode 100644 index 0000000..9cf1864 --- /dev/null +++ b/models/Colbert_HighRes_Smile.STL diff --git a/models/film.stl b/models/film.stl new file mode 100644 index 0000000..30c0f65 --- /dev/null +++ b/models/film.stl @@ -0,0 +1,86 @@ +solid OpenSCAD_Model + facet normal -1.000000 0.000000 0.000000 + outer loop + vertex -0.005000 -1.000000 1.000000 + vertex -0.005000 1.000000 1.000000 + vertex -0.005000 -1.000000 -1.000000 + endloop + endfacet + facet normal -1.000000 -0.000000 0.000000 + outer loop + vertex -0.005000 -1.000000 -1.000000 + vertex -0.005000 1.000000 1.000000 + vertex -0.005000 1.000000 -1.000000 + endloop + endfacet + facet normal 0.000000 0.000000 1.000000 + outer loop + vertex -0.005000 -1.000000 1.000000 + vertex 0.005000 -1.000000 1.000000 + vertex 0.005000 1.000000 1.000000 + endloop + endfacet + facet normal 0.000000 -0.000000 1.000000 + outer loop + vertex -0.005000 1.000000 1.000000 + vertex -0.005000 -1.000000 1.000000 + vertex 0.005000 1.000000 1.000000 + endloop + endfacet + facet normal -0.000000 -1.000000 0.000000 + outer loop + vertex -0.005000 -1.000000 -1.000000 + vertex 0.005000 -1.000000 -1.000000 + vertex 0.005000 -1.000000 1.000000 + endloop + endfacet + facet normal 0.000000 -1.000000 0.000000 + outer loop + vertex -0.005000 -1.000000 1.000000 + vertex -0.005000 -1.000000 -1.000000 + vertex 0.005000 -1.000000 1.000000 + endloop + endfacet + facet normal 0.000000 0.000000 -1.000000 + outer loop + vertex -0.005000 1.000000 -1.000000 + vertex 0.005000 1.000000 -1.000000 + vertex -0.005000 -1.000000 -1.000000 + endloop + endfacet + facet normal -0.000000 0.000000 -1.000000 + outer loop + vertex -0.005000 -1.000000 -1.000000 + vertex 0.005000 1.000000 -1.000000 + vertex 0.005000 -1.000000 -1.000000 + endloop + endfacet + facet normal 0.000000 1.000000 -0.000000 + outer loop + vertex -0.005000 1.000000 1.000000 + vertex 0.005000 1.000000 1.000000 + vertex -0.005000 1.000000 -1.000000 + endloop + endfacet + facet normal 0.000000 1.000000 0.000000 + outer loop + vertex -0.005000 1.000000 -1.000000 + vertex 0.005000 1.000000 1.000000 + vertex 0.005000 1.000000 -1.000000 + endloop + endfacet + facet normal 1.000000 0.000000 0.000000 + outer loop + vertex 0.005000 -1.000000 -1.000000 + vertex 0.005000 1.000000 -1.000000 + vertex 0.005000 1.000000 1.000000 + endloop + endfacet + facet normal 1.000000 0.000000 -0.000000 + outer loop + vertex 0.005000 -1.000000 1.000000 + vertex 0.005000 -1.000000 -1.000000 + vertex 0.005000 1.000000 1.000000 + endloop + endfacet +endsolid OpenSCAD_Model diff --git a/models/liberty.stl b/models/liberty.stl Binary files differnew file mode 100644 index 0000000..89c7e87 --- /dev/null +++ b/models/liberty.stl diff --git a/photon_map.py b/photon_map.py new file mode 100644 index 0000000..cc07db0 --- /dev/null +++ b/photon_map.py @@ -0,0 +1,86 @@ +import numpy as np +from solid import Solid +from geometry import Geometry +from mesh import mesh_from_stl +from sample import uniform_sphere +from pycuda import autoinit +from pycuda.compiler import SourceModule +from pycuda import gpuarray +import pycuda.driver as cuda +import src +import models +from materials import * + +print 'device %s' % autoinit.device.name() + +lens_mesh = mesh_from_stl(models.dir + '/lens.stl') +lens_solid = Solid(0, lens_mesh, material1=pmt_glass, material2=pmt_vacuum) + +film_mesh = mesh_from_stl(models.dir + '/film.stl') +film_solid = Solid(0, film_mesh, displacement=(-5,0,0), material1=h2o, material2=pmt_vacuum) + +sphere_mesh = mesh_from_stl(models.dir + '/sphere.stl') +sphere_solid = Solid(0, sphere_mesh, displacement=(10,0,0), material1=h2o, material2=pmt_vacuum, surface=aluminum) + +geometry = Geometry() +geometry.add_solid(lens_solid) +geometry.add_solid(film_solid) +geometry.add_solid(sphere_solid) +geometry.build(bits=8) + +module = SourceModule(src.kernel, options=['-I' + src.dir], no_extern_c=True, cache_dir=False) +geometry.load(module) +cuda_propagate = module.get_function('propagate') + +nphotons = 100 + +positions = np.tile((10,5,0), nphotons).reshape(nphotons,3) +positions_float3 = np.empty(positions.shape[0], dtype=gpuarray.vec.float3) +positions_float3['x'] = positions[:,0] +positions_float3['y'] = positions[:,1] +positions_float3['z'] = positions[:,2] +positions_gpu = cuda.to_device(positions_float3) + +directions = uniform_sphere(nphotons) +directions_float3 = np.empty(directions.shape[0], dtype=gpuarray.vec.float3) +directions_float3['x'] = directions[:,0] +directions_float3['y'] = directions[:,1] +directions_float3['z'] = directions[:,2] +directions_gpu = cuda.to_device(directions_float3) + +wavelengths = np.random.uniform(200, 800, nphotons).astype(np.float32) +wavelengths_gpu = cuda.to_device(wavelengths) + +times = np.tile(0, nphotons).astype(np.float32) +times_gpu = cuda.to_device(times) + +polarizations = uniform_sphere(nphotons) +polarizations_float3 = np.empty(polarizations.shape[0], dtype=gpuarray.vec.float3) +polarizations_float3['x'] = polarizations[:,0] +polarizations_float3['y'] = polarizations[:,1] +polarizations_float3['z'] = polarizations[:,2] +polarizations_gpu = cuda.to_device(polarizations_float3) + +last_hit_triangles = np.empty(nphotons, dtype=np.int32) +last_hit_triangles_gpu = cuda.to_device(last_hit_triangles) + +states = np.empty(nphotons, dtype=np.int32) +states_gpu = cuda.to_device(states) + +nblocks = 64 + +gpu_kwargs = {'block': (nblocks,1,1), 'grid': (nphotons/nblocks+1,1)} + +cuda_propagate(np.int32(nphotons), positions_gpu, directions_gpu, wavelengths_gpu, polarizations_gpu, times_gpu, states_gpu, last_hit_triangles_gpu, np.int32(geometry.node_map.size-1), np.int32(geometry.first_node), **gpu_kwargs) + +cuda.memcpy_dtoh(states, states_gpu) + +print states + + + + +#from view import view + +#view(geometry) + @@ -1,9 +1,7 @@ import numpy as np class Solid(object): - def __init__(self, id, mesh, rotation=np.identity(3), displacement=(0,0,0), - material1=None, material2=None, surface1=None, surface2=None, - color=0xffffffff): + def __init__(self, id, mesh, rotation=np.identity(3), displacement=(0,0,0), material1=None, material2=None, surface=None, color=0xffffffff): self.id = id self.mesh = mesh @@ -33,19 +31,12 @@ class Solid(object): else: self.material2 = np.tile(material2, len(self.mesh)) - if np.iterable(surface1): - if len(surface1) != len(mesh): + if np.iterable(surface): + if len(surface) != len(mesh): raise ValueError('shape mismatch') - self.surface1 = np.array(surface1, dtype=np.object) + self.surface = np.array(surface, dtype=np.object) else: - self.surface1 = np.tile(surface1, len(self.mesh)) - - if np.iterable(surface2): - if len(surface2) != len(mesh): - raise ValueError('shape mismatch') - self.surface2 = np.array(surface2, dtype=np.object) - else: - self.surface2 = np.tile(surface2, len(self.mesh)) + self.surface = np.tile(surface, len(self.mesh)) if np.iterable(color): if len(color) != len(mesh): diff --git a/src/intersect.h b/src/intersect.h index 92bcf0c..e6566f3 100644 --- a/src/intersect.h +++ b/src/intersect.h @@ -1,4 +1,8 @@ //-*-c-*- + +#ifndef __INTERSECT_H__ +#define __INTERSECT_H__ + #include <math_constants.h> #include "linalg.h" @@ -145,3 +149,5 @@ __device__ bool intersect_box(const float3 &origin, const float3 &direction, con return true; } + +#endif diff --git a/src/kernel.cu b/src/kernel.cu index 796de54..b2e0903 100644 --- a/src/kernel.cu +++ b/src/kernel.cu @@ -7,18 +7,20 @@ #include "rotate.h" #include "intersect.h" #include "materials.h" +#include "physical_constants.h" +#include "random.h" #define STACK_SIZE 500 +#define MAX_DEPTH 5 /* flattened triangle mesh */ texture<float4, 1, cudaReadModeElementType> vertices; __device__ uint4 *triangles; /* material/surface index lookup for each triangle */ -texture<int, 1, cudaReadModeElementType> material1_index; -texture<int, 1, cudaReadModeElementType> material2_index; -texture<int, 1, cudaReadModeElementType> surface1_index; -texture<int, 1, cudaReadModeElementType> surface2_index; +texture<int, 1, cudaReadModeElementType> material1_lookup; +texture<int, 1, cudaReadModeElementType> material2_lookup; +texture<int, 1, cudaReadModeElementType> surface_lookup; /* lower/upper bounds for the bounding box associated with each node/leaf */ texture<float4, 1, cudaReadModeElementType> upper_bounds; @@ -48,11 +50,10 @@ __device__ bool intersect_node(const float3 &origin, const float3 &direction, co direction `direction` and the global mesh texture. If the ray intersects the texture return the index of the triangle which the ray intersected, else return -1. */ -__device__ int intersect_mesh(const float3 &origin, const float3& direction, const int start_node, const int first_node) +__device__ int intersect_mesh(const float3 &origin, const float3& direction, const int start_node, const int first_node, float &distance) { int triangle_idx = -1; - float distance; float min_distance; if (!intersect_node(origin, direction, start_node)) @@ -133,20 +134,8 @@ __global__ void set_pointer(uint4 *triangle_ptr) /* Initialize random number states */ __global__ void init_rng(unsigned long long seed, unsigned long long offset) { - int idx = blockIdx.x*blockDim.x + threadIdx.x; - curand_init(seed, idx, offset, rng_states+idx); -} - -/* */ -__global__ void uniform_sphere(int nthreads, float3 *points) -{ - int idx = blockIdx.x*blockDim.x + threadIdx.x; - - if (idx >= nthreads) - return; - - //curandState rng = *(rng_states+idx); - + int id = blockIdx.x*blockDim.x + threadIdx.x; + curand_init(seed, id, offset, rng_states+id); } /* Translate `points` by the vector `v` */ @@ -188,7 +177,9 @@ __global__ void ray_trace(int nthreads, float3 *origins, float3 *directions, int float3 direction = *(directions+idx); direction /= norm(direction); - int intersection_idx = intersect_mesh(origin, direction, start_node, first_node); + float distance; + + int intersection_idx = intersect_mesh(origin, direction, start_node, first_node, distance); if (intersection_idx == -1) { @@ -206,23 +197,244 @@ __global__ void ray_trace(int nthreads, float3 *origins, float3 *directions, int } } // ray_trace -/* Propagate the photons starting at `origins` traveling in the direction +/* Propagate the photons starting at `positions` traveling in the direction `directions` to their intersection with the global mesh. If the ray intersects the mesh set the hit_solid array value associated with the photon to the triangle index of the triangle the photon intersected, else set the hit_solid array value to -1. */ -__global__ void propagate(int nthreads, float3 *origins, float3 *directions, int start_node, int first_node, int *hit_triangles) +__global__ void propagate(int nthreads, float3 *positions, float3 *directions, float *wavelengths, float3 *polarizations, float *times, int *states, int *last_hit_triangles, int start_node, int first_node) { - int idx = blockIdx.x*blockDim.x + threadIdx.x; + int id = blockIdx.x*blockDim.x + threadIdx.x; - if (idx >= nthreads) + //states[id] = interp_property(materials[0].refractive_index + + //return; + + if (id >= nthreads) return; - float3 origin = *(origins+idx); - float3 direction = *(directions+idx); + float3 position = positions[id]; + float3 direction = directions[id]; direction /= norm(direction); + float3 polarization = polarizations[id]; + float wavelength = wavelengths[id]; + float time = times[id]; + + int last_hit_triangle; + + curandState rng = rng_states[id]; + + unsigned int depth=0; + while (true) + { + depth++; + + if (depth > MAX_DEPTH) + { + states[id] = MAX_DEPTH_REACHED; + break; + } + + float distance; + + last_hit_triangle = intersect_mesh(position, direction, start_node, first_node, distance); + + if (last_hit_triangle == -1) + { + states[id] = NO_HIT; + break; + } + + uint4 triangle_data = triangles[last_hit_triangle]; + + float3 v0 = make_float3(tex1Dfetch(vertices, triangle_data.x)); + float3 v1 = make_float3(tex1Dfetch(vertices, triangle_data.y)); + float3 v2 = make_float3(tex1Dfetch(vertices, triangle_data.z)); + + int material_in_index = 0xFF & (triangle_data.w >> 24); + int material_out_index = 0xFF & (triangle_data.w >> 16); + int surface_index = 0xFF & (triangle_data.w >> 8); + + if (material_in_index < 0 || material_out_index < 0 || surface_index < 0) + { + states[id] = -2; + break; + } + + float3 surface_normal = cross(v1-v0,v2-v0); + surface_normal /= norm(surface_normal); + + Material material1, material2; + if (dot(surface_normal,-direction) > 0.0f) + { + // outside to inside + material1 = materials[material_out_index]; + material2 = materials[material_in_index]; + } + else + { + // inside to outside + material1 = materials[material_in_index]; + material2 = materials[material_out_index]; + surface_normal = -surface_normal; + } + + float refractive_index1 = interp_property(wavelength, material1.refractive_index); + float refractive_index2 = interp_property(wavelength, material2.refractive_index); + + float absorption_length = interp_property(wavelength, material1.absorption_length); + float scattering_length = interp_property(wavelength, material1.scattering_length); + + float absorption_distance = -absorption_length*logf(curand_uniform(&rng)); + float scattering_distance = -absorption_length*logf(curand_uniform(&rng)); + + if (absorption_distance <= scattering_distance) + { + if (absorption_distance <= distance) + { + time = absorption_distance/(SPEED_OF_LIGHT/refractive_index1); + position = position + absorption_distance*direction; + states[id] = BULK_ABSORB; + break; + } // photon is absorbed in material1 + } + else + { + if (scattering_distance <= distance) + { + time = scattering_distance/(SPEED_OF_LIGHT/refractive_index1); + position = position + scattering_distance*direction; + + float x,y; + while (true) + { + y = curand_uniform(&rng); + x = uniform(0, 2*PI); + + if (y < powf(cosf(x),2)) + break; + } + + float phi = uniform(0, 2*PI); + + float3 old_direction = direction; + float3 old_polarization = polarization; + + direction = rotate(old_direction, x, cross(old_polarization, old_direction)); + polarization = rotate(old_polarization, x, cross(old_polarization, old_direction)); + + direction = rotate(direction, phi, old_polarization); + polarization = rotate(polarization, phi, old_polarization); + + direction /= norm(direction); + polarization /= norm(polarization); + + continue; + } // photon is scattered in material1 + } + + if (surface_index != -1) + { + Surface surface = surfaces[surface_index]; + + float absorption = interp_property(wavelength, surface.absorption); + float reflection_diffuse = interp_property(wavelength, surface.reflection_diffuse); + float reflection_specular = interp_property(wavelength, surface.reflection_specular); + + float sum = absorption + reflection_diffuse + reflection_specular; + absorption /= sum; + reflection_diffuse /= sum; + reflection_specular /= sum; + + float p = curand_uniform(&rng); + + if (p < absorption) + { + // absorb + states[id] = SURFACE_ABSORB; + break; + } + else if (p < absorption + reflection_diffuse) + { + // diffusely reflect + while (true) + { + direction = uniform_sphere(&rng); + + if (dot(direction, surface_normal) > 0.0f) + { + // randomize polarization + polarization = cross(uniform_sphere(&rng), direction); + break; + } + } + + continue; + } + else + { + // specularly reflect + direction = rotate(direction, -PI/2.0f, cross(direction, surface_normal)); + + // polarization ? + + continue; + } + } // surface + + // compute reflection/transmission probability + + // p is normal to the plane of incidence + float3 p = cross(direction, surface_normal); + + float normal_coefficient = dot(polarization, p); + float normal_probability = normal_coefficient*normal_coefficient; + + float incident_angle = acosf(dot(surface_normal,-direction)); + float refracted_angle = asinf(sinf(incident_angle)*refractive_index1/refractive_index2); + + if (curand_uniform(&rng) < normal_probability) + { + // photon polarization normal to plane of incidence + float reflection_coefficient = -sinf(incident_angle-refracted_angle)/sinf(incident_angle+refracted_angle); + float reflection_probability = reflection_coefficient*reflection_coefficient; + + polarization = p; + + if (curand_uniform(&rng) < reflection_probability) + { + // reflect off surface + direction = rotate(surface_normal, PI/2.0f, p); + } + else + { + // transmit + direction = rotate(surface_normal, PI-refracted_angle, p); + } + continue; + } + else + { + // photon polarization parallel to surface + float reflection_coefficient = tanf(incident_angle-refracted_angle)/tanf(incident_angle+refracted_angle); + float reflection_probability = reflection_coefficient*reflection_coefficient; + + if (curand_uniform(&rng) < reflection_probability) + { + // reflect off surface + direction = rotate(surface_normal, PI/2.0f, p); + polarization = cross(p, direction); + } + else + { + // transmit + direction = rotate(surface_normal, PI-refracted_angle, p); + polarization = cross(p, direction); + } + continue; + } - *(hit_triangles+idx) = intersect_mesh(origin, direction, start_node, first_node); + } // while(true) } // propagate diff --git a/src/materials.h b/src/materials.h index 05aa121..47b7d22 100644 --- a/src/materials.h +++ b/src/materials.h @@ -1,3 +1,20 @@ +#ifndef __MATERIALS_H__ +#define __MATERIALS_H__ + +__device__ float min_wavelength; +__device__ float max_wavelength; +__device__ float wavelength_step; +__device__ unsigned int wavelength_size; + +enum +{ + INIT = -1, + MAX_DEPTH_REACHED, + NO_HIT, + BULK_ABSORB, + SURFACE_ABSORB +}; + struct Material { float *refractive_index; @@ -12,12 +29,34 @@ struct Surface float *reflection_specular; }; -__device__ struct Material materials[100]; -__device__ struct Surface surfaces[100]; +__device__ Material materials[100]; +__device__ Surface surfaces[100]; + +__device__ float interp_property(const float &x, const float *fp) +{ + if (x < min_wavelength) + return fp[0]; + + if (x > max_wavelength) + return fp[wavelength_size-1]; + + unsigned int jl = (x-min_wavelength)/wavelength_step; + float xl = min_wavelength + jl*wavelength_step; + + return xl + (x-xl)*(fp[jl+1]-fp[jl])/wavelength_step; +} extern "C" { +__global__ void set_wavelength_range(float _min_wavelength, float _max_wavelength, float _wavelength_step, unsigned int _wavelength_size) +{ + min_wavelength = _min_wavelength; + max_wavelength = _max_wavelength; + wavelength_step = _wavelength_step; + wavelength_size = _wavelength_size; +} + __global__ void set_material(int material_index, float *refractive_index, float *absorption_length, float *scattering_length) { materials[material_index].refractive_index = refractive_index; @@ -33,3 +72,5 @@ __global__ void set_surface(int surface_index, float *absorption, float *reflect } } // extern "c" + +#endif diff --git a/src/physical_constants.h b/src/physical_constants.h new file mode 100644 index 0000000..2ff87cd --- /dev/null +++ b/src/physical_constants.h @@ -0,0 +1,7 @@ +#ifndef __PHYSICAL_CONSTANTS_H__ +#define __PHYSICAL_CONSTANTS_H__ + +#define SPEED_OF_LIGHT 2.99792458e8 +#define PI 3.141592653589793 + +#endif diff --git a/src/random.h b/src/random.h new file mode 100644 index 0000000..74329da --- /dev/null +++ b/src/random.h @@ -0,0 +1,22 @@ +#ifndef __RANDOM_H__ +#define __RANDOM_H__ + +#include <curand_kernel.h> + +#include "physical_constants.h" + +__device__ float uniform(curandState *rng, float a=0.0, float b=1.0) +{ + return a + curand_uniform(rng)*(b-a); +} + +__device__ float3 uniform_sphere(curandState *rng) +{ + float theta = uniform(rng, 0, 2*PI); + float u = uniform(rng, -1, 1); + float c = sqrtf(1-u*u); + + return make_float3(c*cosf(theta), c*sinf(theta), u); +} + +#endif |