summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--camera.py261
-rw-r--r--color/Illuminanta.csv531
-rw-r--r--color/chromaticity.py18
-rw-r--r--color/scvle_1.csv401
-rw-r--r--color/vl1924e_1.csv471
-rw-r--r--detectors/__init__.py4
-rw-r--r--geometry.py14
-rw-r--r--solids/__init__.py6
-rw-r--r--src/kernel.cu250
-rw-r--r--src/photon.h8
10 files changed, 1697 insertions, 267 deletions
diff --git a/camera.py b/camera.py
index 2feafd7..18203c0 100644
--- a/camera.py
+++ b/camera.py
@@ -2,8 +2,10 @@ import numpy as np
from itertools import product, count
from threading import Thread, Lock
import time
+import datetime
import os
import sys
+from color import map_wavelength
from transform import rotate
@@ -14,6 +16,32 @@ from pycuda import gpuarray
from pycuda.characterize import sizeof
import pycuda.driver as cuda
+from subprocess import call
+import shutil
+import tempfile
+
+def timeit(func):
+ def f(*args, **kwargs):
+ t0 = time.time()
+ func(*args, **kwargs)
+ elapsed = time.time() - t0
+ print '%s elapsed in %s().' % (datetime.timedelta(seconds=elapsed), func.__name__)
+ return f
+
+def encode_movie(dir):
+ root, ext = 'movie', 'avi'
+ for i in count():
+ filename = '.'.join([root + str(i).zfill(5), ext])
+
+ if not os.path.exists(filename):
+ break
+
+ call(['mencoder', 'mf://' + dir + '/*.png', '-mf', 'fps=10', '-o', filename, '-ovc', 'xvid', '-xvidencopts', 'bitrate=3000'])
+
+ shutil.rmtree(dir)
+
+ print 'movie saved to %s.' % filename
+
def get_rays(position, size = (800, 600), film_size = (0.035, 0.024), focal_length=0.05):
"""
Generate ray positions and directions from a pinhole camera facing the negative y direction.
@@ -62,9 +90,11 @@ class Camera(Thread):
self.context.push()
self.ray_trace_kernel = self.module.get_function('ray_trace')
self.rotate_kernel = self.module.get_function('rotate')
+ self.rotate_around_point_kernel = self.module.get_function('rotate_around_point')
self.translate_kernel = self.module.get_function('translate')
- self.build_rgb_lookup_kernel = self.module.get_function('build_rgb_lookup')
- self.render_kernel = self.module.get_function('render')
+ self.update_xyz_lookup_kernel = self.module.get_function('update_xyz_lookup')
+ self.update_xyz_image_kernel = self.module.get_function('update_xyz_image')
+ self.process_image_kernel = self.module.get_function('process_image')
self.init_rng_kernel = self.module.get_function('init_rng')
self.context.pop()
@@ -77,7 +107,7 @@ class Camera(Thread):
self.axis1 = np.array([1,0,0], float)
self.axis2 = np.array([0,0,1], float)
- origins, directions = get_rays(self.point)
+ origins, directions = get_rays(self.point, self.size)
with self.lock:
self.context.push()
@@ -86,85 +116,84 @@ class Camera(Thread):
self.pixels_gpu = gpuarray.zeros(self.width*self.height, dtype=np.int32)
self.context.pop()
+ self.movie = False
+ self.movie_index = 0
+ self.movie_dir = None
self.render = False
- def build_rgb_lookup(self, source_position=(0,0,0)):
+ @timeit
+ def initialize_render(self):
with self.lock:
self.context.push()
- print '\ninitializing random number states.'
self.rng_states_gpu = cuda.mem_alloc(self.width*self.height*sizeof('curandStateXORWOW', '#include <curand_kernel.h>'))
self.init_rng_kernel(np.int32(self.width*self.height), self.rng_states_gpu, np.int32(0), np.int32(0), block=(self.nblocks,1,1), grid=(self.width*self.height//self.nblocks+1,1))
+ self.xyz_lookup1_gpu = gpuarray.zeros(len(self.geometry.mesh.triangles), dtype=gpuarray.vec.float3)
+ self.xyz_lookup2_gpu = gpuarray.zeros(len(self.geometry.mesh.triangles), dtype=gpuarray.vec.float3)
+ self.image_gpu = gpuarray.zeros(self.width*self.height, dtype=gpuarray.vec.float3)
self.context.synchronize()
- print 'done.'
-
- self.source_positions_gpu = gpuarray.empty(len(geometry.mesh.triangles), dtype=gpuarray.vec.float3)
- self.source_positions_gpu.fill(gpuarray.vec.make_float3(*source_position))
+ self.context.pop()
- source_directions = np.mean(self.geometry.mesh[:], axis=1) - source_position
- self.source_directions_gpu = gpuarray.to_gpu(source_directions.astype(np.float32).view(gpuarray.vec.float3))
+ self.source_position = self.point
- self.rgb_lookup1_gpu = gpuarray.zeros(self.source_positions_gpu.size, dtype=gpuarray.vec.float3)
- self.rgb_lookup2_gpu = gpuarray.zeros(self.source_positions_gpu.size, dtype=gpuarray.vec.float3)
+ self.nimages = 0
+ self.nlookup_calls = 0
+ self.max_steps = 10
- self.max_steps = 10
- rgb_runs = 100
-
- print 'building rgb lookup.'
- self.build_rgb_lookup_kernel(np.int32(self.source_positions_gpu.size), self.source_positions_gpu, self.source_directions_gpu, self.rgb_lookup1_gpu, self.rgb_lookup2_gpu, np.uint64(np.random.randint(np.iinfo(np.int64).max)), np.int32(rgb_runs), np.int32(self.max_steps), block=(self.nblocks,1,1), grid=(self.source_positions_gpu.size//self.nblocks+1,1))
- self.context.synchronize()
- print 'done.'
+ def clear_xyz_lookup(self):
+ with self.lock:
+ self.context.push()
+ self.xyz_lookup1_gpu.fill(gpuarray.vec.make_float3(0.0,0.0,0.0))
+ self.xyz_lookup2_gpu.fill(gpuarray.vec.make_float3(0.0,0.0,0.0))
+ self.context.pop()
- rgb_lookup1 = self.rgb_lookup1_gpu.get().view(np.float32)
- rgb_lookup1 /= rgb_runs
- rgb_lookup1[rgb_lookup1 > 1.0] = 1.0
- self.rgb_lookup1_gpu.set(rgb_lookup1.view(gpuarray.vec.float3))
+ self.nlookup_calls = 0
- rgb_lookup2 = self.rgb_lookup2_gpu.get().view(np.float32)
- rgb_lookup2 /= rgb_runs
- rgb_lookup2[rgb_lookup2 > 1.0] = 1.0
- self.rgb_lookup2_gpu.set(rgb_lookup2.view(gpuarray.vec.float3))
+ def update_xyz_lookup(self, source_position):
+ with self.lock:
+ self.context.push()
+ for i in range(self.xyz_lookup1_gpu.size//(self.width*self.height)+1):
+ self.update_xyz_lookup_kernel(np.int32(self.width*self.height), np.int32(self.xyz_lookup1_gpu.size), np.int32(i*self.width*self.height), gpuarray.vec.make_float3(*source_position), self.rng_states_gpu, np.float32(685.0), gpuarray.vec.make_float3(1.0,0.0,0.0), self.xyz_lookup1_gpu, self.xyz_lookup2_gpu, np.int32(self.max_steps), block=(self.nblocks,1,1), grid=(self.width*self.height//self.nblocks+1,1))
- self.nimages = 0
- self.rgb_image = np.zeros((self.pixels_gpu.size, 3), float)
- self.r_gpu = gpuarray.zeros(self.pixels_gpu.size, dtype=np.float32)
- self.g_gpu = gpuarray.zeros(self.pixels_gpu.size, dtype=np.float32)
- self.b_gpu = gpuarray.zeros(self.pixels_gpu.size, dtype=np.float32)
+ for i in range(self.xyz_lookup1_gpu.size//(self.width*self.height)+1):
+ self.update_xyz_lookup_kernel(np.int32(self.width*self.height), np.int32(self.xyz_lookup1_gpu.size), np.int32(i*self.width*self.height), gpuarray.vec.make_float3(*source_position), self.rng_states_gpu, np.float32(545.0), gpuarray.vec.make_float3(0.0,1.0,0.0), self.xyz_lookup1_gpu, self.xyz_lookup2_gpu, np.int32(self.max_steps), block=(self.nblocks,1,1), grid=(self.width*self.height//self.nblocks+1,1))
+ for i in range(self.xyz_lookup1_gpu.size//(self.width*self.height)+1):
+ self.update_xyz_lookup_kernel(np.int32(self.width*self.height), np.int32(self.xyz_lookup1_gpu.size), np.int32(i*self.width*self.height), gpuarray.vec.make_float3(*source_position), self.rng_states_gpu, np.float32(445.0), gpuarray.vec.make_float3(0.0,0.0,1.0), self.xyz_lookup1_gpu, self.xyz_lookup2_gpu, np.int32(self.max_steps), block=(self.nblocks,1,1), grid=(self.width*self.height//self.nblocks+1,1))
self.context.pop()
+ self.nlookup_calls += 1
+
def clear_image(self):
- try:
- self.rgb_image.fill(0.0)
- self.nimages = 0
- except AttributeError:
- pass
+ with self.lock:
+ self.context.push()
+ self.image_gpu.fill(gpuarray.vec.make_float3(0.0,0.0,0.0))
+ self.context.pop()
- def acquire_image(self):
- if not hasattr(self, 'rng_states_gpu'):
- self.build_rgb_lookup()
+ self.nimages = 0
+ def update_image(self):
with self.lock:
self.context.push()
- self.render_kernel(np.int32(self.pixels_gpu.size), self.rng_states_gpu, self.origins_gpu, self.directions_gpu, self.r_gpu, self.g_gpu, self.b_gpu, self.rgb_lookup1_gpu, self.rgb_lookup2_gpu, np.int32(self.max_steps), block=(self.nblocks,1,1), grid=(self.pixels_gpu.size//self.nblocks+1,1))
- self.rgb_image[:,0] += self.r_gpu.get()
- self.rgb_image[:,1] += self.g_gpu.get()
- self.rgb_image[:,2] += self.b_gpu.get()
- self.nimages += 1
- self.context.pop()
+ self.update_xyz_image_kernel(np.int32(self.width*self.height), self.rng_states_gpu, self.origins_gpu, self.directions_gpu, np.float32(685.0), gpuarray.vec.make_float3(1.0,0.0,0.0), self.xyz_lookup1_gpu, self.xyz_lookup2_gpu, self.image_gpu, np.int32(self.nlookup_calls), np.int32(self.max_steps), block=(self.nblocks,1,1), grid=(self.width*self.height//self.nblocks+1,1))
+
+ self.update_xyz_image_kernel(np.int32(self.width*self.height), self.rng_states_gpu, self.origins_gpu, self.directions_gpu, np.float32(545.0), gpuarray.vec.make_float3(0.0,1.0,0.0), self.xyz_lookup1_gpu, self.xyz_lookup2_gpu, self.image_gpu, np.int32(self.nlookup_calls), np.int32(self.max_steps), block=(self.nblocks,1,1), grid=(self.width*self.height//self.nblocks+1,1))
- def render_image(self):
- scaled_rgb_image = ((self.rgb_image/self.nimages)*255).astype(np.int32)
+ self.update_xyz_image_kernel(np.int32(self.width*self.height), self.rng_states_gpu, self.origins_gpu, self.directions_gpu, np.float32(445.0), gpuarray.vec.make_float3(0.0,0.0,1.0), self.xyz_lookup1_gpu, self.xyz_lookup2_gpu, self.image_gpu, np.int32(self.nlookup_calls), np.int32(self.max_steps), block=(self.nblocks,1,1), grid=(self.width*self.height//self.nblocks+1,1))
+ self.context.pop()
- image = scaled_rgb_image[:,0] << 16 | scaled_rgb_image[:,1] << 8 | scaled_rgb_image[:,2]
+ self.nimages += 1
- pygame.surfarray.blit_array(self.screen, image.reshape(self.size))
- pygame.display.flip()
+ def process_image(self):
+ with self.lock:
+ self.context.push()
+ self.process_image_kernel(np.int32(self.width*self.height), self.image_gpu, self.pixels_gpu, np.int32(self.nimages), block=(self.nblocks,1,1), grid=((self.width*self.height)//self.nblocks+1,1))
+ self.context.pop()
- def screenshot(self, dir=''):
+ def screenshot(self, dir='', start=0):
root, ext = 'screenshot', 'png'
- for i in count():
- filename = os.path.join(dir, '.'.join([root + str(i).zfill(4), ext]))
+ for i in count(start):
+ filename = os.path.join(dir, '.'.join([root + str(i).zfill(5), ext]))
if not os.path.exists(filename):
break
@@ -175,14 +204,34 @@ class Camera(Thread):
def rotate(self, phi, n):
with self.lock:
self.context.push()
- self.rotate_kernel(np.int32(self.pixels_gpu.size), self.origins_gpu, phi, gpuarray.vec.make_float3(*n), block=(self.nblocks,1,1), grid=(self.pixels_gpu.size//self.nblocks+1,1))
- self.rotate_kernel(np.int32(self.pixels_gpu.size), self.directions_gpu, phi, gpuarray.vec.make_float3(*n), block=(self.nblocks,1,1), grid=(self.pixels_gpu.size//self.nblocks+1,1))
+ self.rotate_kernel(np.int32(self.pixels_gpu.size), self.origins_gpu, np.float32(phi), gpuarray.vec.make_float3(*n), block=(self.nblocks,1,1), grid=(self.pixels_gpu.size//self.nblocks+1,1))
+ self.rotate_kernel(np.int32(self.pixels_gpu.size), self.directions_gpu, np.float32(phi), gpuarray.vec.make_float3(*n), block=(self.nblocks,1,1), grid=(self.pixels_gpu.size//self.nblocks+1,1))
self.point = rotate(self.point, phi, n)
self.axis1 = rotate(self.axis1, phi, n)
self.axis2 = rotate(self.axis2, phi, n)
self.context.pop()
+ if self.render:
+ self.clear_image()
+
+ self.update()
+
+ def rotate_around_point(self, phi, n, point):
+ with self.lock:
+ self.context.push()
+ self.rotate_around_point_kernel(np.int32(self.origins_gpu.size), self.origins_gpu, np.float32(phi), gpuarray.vec.make_float3(*n), gpuarray.vec.make_float3(*point), block=(self.nblocks,1,1), grid=(self.origins_gpu.size//self.nblocks+1,1))
+ self.rotate_kernel(np.int32(self.directions_gpu.size), self.directions_gpu, np.float32(phi), gpuarray.vec.make_float3(*n), block=(self.nblocks,1,1), grid=(self.directions_gpu.size//self.nblocks+1,1))
+ self.context.pop()
+
+ self.axis1 = rotate(self.axis1, phi, n)
+ self.axis2 = rotate(self.axis2, phi, n)
+
+ if self.render:
+ self.clear_image()
+
+ self.update()
+
def translate(self, v):
with self.lock:
self.context.push()
@@ -191,55 +240,58 @@ class Camera(Thread):
self.point += v
self.context.pop()
+ if self.render:
+ self.clear_image()
+
+ self.update()
+
def update(self):
if self.render:
- self.acquire_image()
- self.render_image()
- return
+ while self.nlookup_calls < 10:
+ self.update_xyz_lookup(self.source_position)
+ self.update_image()
+ self.process_image()
+ else:
+ with self.lock:
+ self.context.push()
+ self.ray_trace_kernel(np.int32(self.pixels_gpu.size), self.origins_gpu, self.directions_gpu, self.pixels_gpu, block=(self.nblocks,1,1), grid=(self.pixels_gpu.size//self.nblocks+1,1))
+ self.context.pop()
with self.lock:
self.context.push()
- t0 = time.time()
- self.ray_trace_kernel(np.int32(self.pixels_gpu.size), self.origins_gpu, self.directions_gpu, self.pixels_gpu, block=(self.nblocks,1,1), grid=(self.pixels_gpu.size//self.nblocks+1,1))
- self.context.synchronize()
- elapsed = time.time() - t0
-
- print '\relapsed %f sec.' % elapsed,
- sys.stdout.flush()
-
pygame.surfarray.blit_array(self.screen, self.pixels_gpu.get().reshape(self.size))
pygame.display.flip()
self.context.pop()
+ if self.movie:
+ self.screenshot(self.movie_dir, self.movie_index)
+ self.movie_index += 1
+
def run(self):
self.update()
done = False
clicked = False
shift = False
+ ctrl = False
- current_layer = 0
+ #current_layer = 0
while not done:
self.clock.tick(20)
if self.render and not clicked and not pygame.event.peek(KEYDOWN):
- self.acquire_image()
- self.render_image()
+ self.update()
for event in pygame.event.get():
if event.type == MOUSEBUTTONDOWN:
if event.button == 4:
v = self.scale*np.cross(self.axis1,self.axis2)/10.0
self.translate(v)
- self.clear_image()
- self.update()
if event.button == 5:
v = -self.scale*np.cross(self.axis1,self.axis2)/10.0
self.translate(v)
- self.clear_image()
- self.update()
if event.button == 1:
clicked = True
@@ -263,56 +315,56 @@ class Camera(Thread):
if shift:
v = mouse_direction*self.scale*length/float(self.width)
self.translate(v)
- self.clear_image()
- self.update()
else:
phi = np.float32(2*np.pi*length/float(self.width))
n = rotate(mouse_direction, np.pi/2, -np.cross(self.axis1,self.axis2))
- self.rotate(phi, n)
- self.clear_image()
- self.update()
+ if ctrl:
+ self.rotate_around_point(phi, n, self.point)
+ else:
+ self.rotate(phi, n)
if event.type == KEYDOWN:
if event.key == K_a:
v = self.scale*self.axis1/10.0
self.translate(v)
- self.clear_image()
- self.update()
if event.key == K_d:
v = -self.scale*self.axis1/10.0
self.translate(v)
- self.clear_image()
- self.update()
if event.key == K_w:
v = self.scale*np.cross(self.axis1,self.axis2)/10.0
self.translate(v)
- self.clear_image()
- self.update()
if event.key == K_s:
v = -self.scale*np.cross(self.axis1,self.axis2)/10.0
self.translate(v)
- self.clear_image()
- self.update()
if event.key == K_SPACE:
v = self.scale*self.axis2/10.0
self.translate(v)
- self.clear_image()
- self.update()
- if event.key == K_LCTRL:
- v = -self.scale*self.axis2/10.0
- self.translate(v)
+ # if event.key == K_LCTRL:
+ # v = -self.scale*self.axis2/10.0
+ # self.translate(v)
+
+ if event.key == K_F6:
+ self.clear_xyz_lookup()
self.clear_image()
- self.update()
+ self.source_position = self.point
+
+ if event.key == K_p:
+ for i in range(10):
+ self.update_xyz_lookup(self.point)
+ self.source_position = self.point
if event.key == K_LSHIFT or event.key == K_RSHIFT:
shift = True
+ if event.key == K_LCTRL or event.key == K_RCTRL:
+ ctrl = True
+
if event.key == K_ESCAPE:
done = True
break
@@ -347,18 +399,37 @@ class Camera(Thread):
self.screenshot()
if event.key == K_F5:
+ if not hasattr(self, 'rng_states_gpu'):
+ self.initialize_render()
+
self.render = not self.render
+ self.clear_image()
self.update()
+ if event.key == K_m:
+ if self.movie:
+ encode_movie(self.movie_dir)
+ self.movie_dir = None
+ self.movie = False
+ else:
+ self.movie_index = 0
+ self.movie_dir = tempfile.mkdtemp()
+ self.movie = True
+
if event.type == KEYUP:
if event.key == K_LSHIFT or event.key == K_RSHIFT:
shift = False
+ if event.key == K_LCTRL or event.key == K_RCTRL:
+ ctrl = False
+
if event.type == pygame.QUIT:
done = True
break
- print
+ if self.movie:
+ encode_movie(self.movie_dir)
+
pygame.display.quit()
if __name__ == '__main__':
diff --git a/color/Illuminanta.csv b/color/Illuminanta.csv
new file mode 100644
index 0000000..3dba86c
--- /dev/null
+++ b/color/Illuminanta.csv
@@ -0,0 +1,531 @@
+300,0.930483
+301,0.967643
+302,1.005970
+303,1.045490
+304,1.086230
+305,1.128210
+306,1.171470
+307,1.216020
+308,1.261880
+309,1.309100
+310,1.357690
+311,1.407680
+312,1.459100
+313,1.511980
+314,1.566330
+315,1.622190
+316,1.679590
+317,1.738550
+318,1.799100
+319,1.861270
+320,1.925080
+321,1.990570
+322,2.057760
+323,2.126670
+324,2.197340
+325,2.269800
+326,2.344060
+327,2.420170
+328,2.498140
+329,2.578010
+330,2.659810
+331,2.743550
+332,2.829280
+333,2.917010
+334,3.006780
+335,3.098610
+336,3.192530
+337,3.288570
+338,3.386760
+339,3.487120
+340,3.589680
+341,3.694470
+342,3.801520
+343,3.910850
+344,4.022500
+345,4.136480
+346,4.252820
+347,4.371560
+348,4.492720
+349,4.616310
+350,4.742380
+351,4.870950
+352,5.002040
+353,5.135680
+354,5.271890
+355,5.410700
+356,5.552130
+357,5.696220
+358,5.842980
+359,5.992440
+360,6.144620
+361,6.299550
+362,6.457240
+363,6.617740
+364,6.781050
+365,6.947200
+366,7.116210
+367,7.288110
+368,7.462920
+369,7.640660
+370,7.821350
+371,8.005010
+372,8.191670
+373,8.381340
+374,8.574040
+375,8.769800
+376,8.968640
+377,9.170560
+378,9.375610
+379,9.583780
+380,9.795100
+381,10.009600
+382,10.227300
+383,10.448100
+384,10.672200
+385,10.899600
+386,11.130200
+387,11.364000
+388,11.601200
+389,11.841600
+390,12.085300
+391,12.332400
+392,12.582800
+393,12.836600
+394,13.093800
+395,13.354300
+396,13.618200
+397,13.885500
+398,14.156300
+399,14.430400
+400,14.708000
+401,14.989100
+402,15.273600
+403,15.561600
+404,15.853000
+405,16.148000
+406,16.446400
+407,16.748400
+408,17.053800
+409,17.362800
+410,17.675300
+411,17.991300
+412,18.310800
+413,18.633900
+414,18.960500
+415,19.290700
+416,19.624400
+417,19.961700
+418,20.302600
+419,20.647000
+420,20.995000
+421,21.346500
+422,21.701600
+423,22.060300
+424,22.422500
+425,22.788300
+426,23.157700
+427,23.530700
+428,23.907200
+429,24.287300
+430,24.670900
+431,25.058100
+432,25.448900
+433,25.843200
+434,26.241100
+435,26.642500
+436,27.047500
+437,27.456000
+438,27.868100
+439,28.283600
+440,28.702700
+441,29.125300
+442,29.551500
+443,29.981100
+444,30.414200
+445,30.850800
+446,31.290900
+447,31.734500
+448,32.181500
+449,32.632000
+450,33.085900
+451,33.543200
+452,34.004000
+453,34.468200
+454,34.935800
+455,35.406800
+456,35.881100
+457,36.358800
+458,36.839900
+459,37.324300
+460,37.812100
+461,38.303100
+462,38.797500
+463,39.295100
+464,39.796000
+465,40.300200
+466,40.807600
+467,41.318200
+468,41.832000
+469,42.349100
+470,42.869300
+471,43.392600
+472,43.919200
+473,44.448800
+474,44.981600
+475,45.517400
+476,46.056300
+477,46.598300
+478,47.143300
+479,47.691300
+480,48.242300
+481,48.796300
+482,49.353300
+483,49.913200
+484,50.476000
+485,51.041800
+486,51.610400
+487,52.181800
+488,52.756100
+489,53.333200
+490,53.913200
+491,54.495800
+492,55.081300
+493,55.669400
+494,56.260300
+495,56.853900
+496,57.450100
+497,58.048900
+498,58.650400
+499,59.254500
+500,59.861100
+501,60.470300
+502,61.082000
+503,61.696200
+504,62.312800
+505,62.932000
+506,63.553500
+507,64.177500
+508,64.803800
+509,65.432500
+510,66.063500
+511,66.696800
+512,67.332400
+513,67.970200
+514,68.610200
+515,69.252500
+516,69.896900
+517,70.543500
+518,71.192200
+519,71.843000
+520,72.495900
+521,73.150800
+522,73.807700
+523,74.466600
+524,75.127500
+525,75.790300
+526,76.455100
+527,77.121700
+528,77.790200
+529,78.460500
+530,79.132600
+531,79.806500
+532,80.482100
+533,81.159500
+534,81.838600
+535,82.519300
+536,83.201700
+537,83.885600
+538,84.571200
+539,85.258400
+540,85.947000
+541,86.637200
+542,87.328800
+543,88.021900
+544,88.716500
+545,89.412400
+546,90.109700
+547,90.808300
+548,91.508200
+549,92.209500
+550,92.912000
+551,93.615700
+552,94.320600
+553,95.026700
+554,95.733900
+555,96.442300
+556,97.151800
+557,97.862300
+558,98.573900
+559,99.286400
+560,100.000000
+561,100.715000
+562,101.430000
+563,102.146000
+564,102.864000
+565,103.582000
+566,104.301000
+567,105.020000
+568,105.741000
+569,106.462000
+570,107.184000
+571,107.906000
+572,108.630000
+573,109.354000
+574,110.078000
+575,110.803000
+576,111.529000
+577,112.255000
+578,112.982000
+579,113.709000
+580,114.436000
+581,115.164000
+582,115.893000
+583,116.622000
+584,117.351000
+585,118.080000
+586,118.810000
+587,119.540000
+588,120.270000
+589,121.001000
+590,121.731000
+591,122.462000
+592,123.193000
+593,123.924000
+594,124.655000
+595,125.386000
+596,126.118000
+597,126.849000
+598,127.580000
+599,128.312000
+600,129.043000
+601,129.774000
+602,130.505000
+603,131.236000
+604,131.966000
+605,132.697000
+606,133.427000
+607,134.157000
+608,134.887000
+609,135.617000
+610,136.346000
+611,137.075000
+612,137.804000
+613,138.532000
+614,139.260000
+615,139.988000
+616,140.715000
+617,141.441000
+618,142.167000
+619,142.893000
+620,143.618000
+621,144.343000
+622,145.067000
+623,145.790000
+624,146.513000
+625,147.235000
+626,147.957000
+627,148.678000
+628,149.398000
+629,150.117000
+630,150.836000
+631,151.554000
+632,152.271000
+633,152.988000
+634,153.704000
+635,154.418000
+636,155.132000
+637,155.845000
+638,156.558000
+639,157.269000
+640,157.979000
+641,158.689000
+642,159.397000
+643,160.104000
+644,160.811000
+645,161.516000
+646,162.221000
+647,162.924000
+648,163.626000
+649,164.327000
+650,165.028000
+651,165.726000
+652,166.424000
+653,167.121000
+654,167.816000
+655,168.510000
+656,169.203000
+657,169.895000
+658,170.586000
+659,171.275000
+660,171.963000
+661,172.650000
+662,173.335000
+663,174.019000
+664,174.702000
+665,175.383000
+666,176.063000
+667,176.741000
+668,177.419000
+669,178.094000
+670,178.769000
+671,179.441000
+672,180.113000
+673,180.783000
+674,181.451000
+675,182.118000
+676,182.783000
+677,183.447000
+678,184.109000
+679,184.770000
+680,185.429000
+681,186.087000
+682,186.743000
+683,187.397000
+684,188.050000
+685,188.701000
+686,189.350000
+687,189.998000
+688,190.644000
+689,191.288000
+690,191.931000
+691,192.572000
+692,193.211000
+693,193.849000
+694,194.484000
+695,195.118000
+696,195.750000
+697,196.381000
+698,197.009000
+699,197.636000
+700,198.261000
+701,198.884000
+702,199.506000
+703,200.125000
+704,200.743000
+705,201.359000
+706,201.972000
+707,202.584000
+708,203.195000
+709,203.803000
+710,204.409000
+711,205.013000
+712,205.616000
+713,206.216000
+714,206.815000
+715,207.411000
+716,208.006000
+717,208.599000
+718,209.189000
+719,209.778000
+720,210.365000
+721,210.949000
+722,211.532000
+723,212.112000
+724,212.691000
+725,213.268000
+726,213.842000
+727,214.415000
+728,214.985000
+729,215.553000
+730,216.120000
+731,216.684000
+732,217.246000
+733,217.806000
+734,218.364000
+735,218.920000
+736,219.473000
+737,220.025000
+738,220.574000
+739,221.122000
+740,221.667000
+741,222.210000
+742,222.751000
+743,223.290000
+744,223.826000
+745,224.361000
+746,224.893000
+747,225.423000
+748,225.951000
+749,226.477000
+750,227.000000
+751,227.522000
+752,228.041000
+753,228.558000
+754,229.073000
+755,229.585000
+756,230.096000
+757,230.604000
+758,231.110000
+759,231.614000
+760,232.115000
+761,232.615000
+762,233.112000
+763,233.606000
+764,234.099000
+765,234.589000
+766,235.078000
+767,235.564000
+768,236.047000
+769,236.529000
+770,237.008000
+771,237.485000
+772,237.959000
+773,238.432000
+774,238.902000
+775,239.370000
+776,239.836000
+777,240.299000
+778,240.760000
+779,241.219000
+780,241.675000
+781,242.130000
+782,242.582000
+783,243.031000
+784,243.479000
+785,243.924000
+786,244.367000
+787,244.808000
+788,245.246000
+789,245.682000
+790,246.116000
+791,246.548000
+792,246.977000
+793,247.404000
+794,247.829000
+795,248.251000
+796,248.671000
+797,249.089000
+798,249.505000
+799,249.918000
+800,250.329000
+801,250.738000
+802,251.144000
+803,251.548000
+804,251.950000
+805,252.350000
+806,252.747000
+807,253.142000
+808,253.535000
+809,253.925000
+810,254.314000
+811,254.700000
+812,255.083000
+813,255.465000
+814,255.844000
+815,256.221000
+816,256.595000
+817,256.968000
+818,257.338000
+819,257.706000
+820,258.071000
+821,258.434000
+822,258.795000
+823,259.154000
+824,259.511000
+825,259.865000
+826,260.217000
+827,260.567000
+828,260.914000
+829,261.259000
+830,261.602000 \ No newline at end of file
diff --git a/color/chromaticity.py b/color/chromaticity.py
index fbd812a..859cc14 100644
--- a/color/chromaticity.py
+++ b/color/chromaticity.py
@@ -1,24 +1,14 @@
import numpy as np
import os
+import sys
dir = os.path.split(os.path.realpath(__file__))[0]
-f = open(dir + '/sbrgb10w.csv')
+sys.path.append(dir + '/..')
-color_map = []
-for line in f:
- color_map.append([float(s) for s in line.split(',')])
+from tools import read_csv
-f.close()
-
-color_map = np.array(color_map)
-
-# zero negative coefficients
-color_map[color_map < 0] = 0
-
-# normalize coefficients
-for i in range(len(color_map)):
- color_map[i,1:] /= np.sum(color_map[i,1:])
+color_map = read_csv(dir + '/ciexyz64_1.csv')
def map_wavelength(wavelength):
r = np.interp(wavelength, color_map[:,0], color_map[:,1])
diff --git a/color/scvle_1.csv b/color/scvle_1.csv
new file mode 100644
index 0000000..0ab7def
--- /dev/null
+++ b/color/scvle_1.csv
@@ -0,0 +1,401 @@
+380, 0.0005890000
+381, 0.0006650000
+382, 0.0007520000
+383, 0.0008540000
+384, 0.0009720000
+385, 0.0011080000
+386, 0.0012680000
+387, 0.0014530000
+388, 0.0016680000
+389, 0.0019180000
+390, 0.0022090000
+391, 0.0025470000
+392, 0.0029390000
+393, 0.0033940000
+394, 0.0039210000
+395, 0.0045300000
+396, 0.0052400000
+397, 0.0060500000
+398, 0.0069800000
+399, 0.0080600000
+400, 0.0092900000
+401, 0.0107000000
+402, 0.0123100000
+403, 0.0141300000
+404, 0.0161900000
+405, 0.0185200000
+406, 0.0211300000
+407, 0.0240500000
+408, 0.0273000000
+409, 0.0308900000
+410, 0.0348400000
+411, 0.0391600000
+412, 0.0439000000
+413, 0.0490000000
+414, 0.0545000000
+415, 0.0604000000
+416, 0.0668000000
+417, 0.0736000000
+418, 0.0808000000
+419, 0.0885000000
+420, 0.0966000000
+421, 0.1052000000
+422, 0.1141000000
+423, 0.1235000000
+424, 0.1334000000
+425, 0.1436000000
+426, 0.1541000000
+427, 0.1651000000
+428, 0.1764000000
+429, 0.1879000000
+430, 0.1998000000
+431, 0.2119000000
+432, 0.2243000000
+433, 0.2369000000
+434, 0.2496000000
+435, 0.2625000000
+436, 0.2755000000
+437, 0.2886000000
+438, 0.3017000000
+439, 0.3149000000
+440, 0.3281000000
+441, 0.3412000000
+442, 0.3543000000
+443, 0.3673000000
+444, 0.3803000000
+445, 0.3931000000
+446, 0.4060000000
+447, 0.4180000000
+448, 0.4310000000
+449, 0.4430000000
+450, 0.4550000000
+451, 0.4670000000
+452, 0.4790000000
+453, 0.4900000000
+454, 0.5020000000
+455, 0.5130000000
+456, 0.5240000000
+457, 0.5350000000
+458, 0.5460000000
+459, 0.5570000000
+460, 0.5670000000
+461, 0.5780000000
+462, 0.5880000000
+463, 0.5990000000
+464, 0.6100000000
+465, 0.6200000000
+466, 0.6310000000
+467, 0.6420000000
+468, 0.6530000000
+469, 0.6640000000
+470, 0.6760000000
+471, 0.6870000000
+472, 0.6990000000
+473, 0.7100000000
+474, 0.7220000000
+475, 0.7340000000
+476, 0.7450000000
+477, 0.7570000000
+478, 0.7690000000
+479, 0.7810000000
+480, 0.7930000000
+481, 0.8050000000
+482, 0.8170000000
+483, 0.8280000000
+484, 0.8400000000
+485, 0.8510000000
+486, 0.8620000000
+487, 0.8730000000
+488, 0.8840000000
+489, 0.8940000000
+490, 0.9040000000
+491, 0.9140000000
+492, 0.9230000000
+493, 0.9320000000
+494, 0.9410000000
+495, 0.9490000000
+496, 0.9570000000
+497, 0.9640000000
+498, 0.9700000000
+499, 0.9760000000
+500, 0.9820000000
+501, 0.9860000000
+502, 0.9900000000
+503, 0.9940000000
+504, 0.9970000000
+505, 0.9980000000
+506, 1.0000000000
+507, 1.0000000000
+508, 1.0000000000
+509, 0.9980000000
+510, 0.9970000000
+511, 0.9940000000
+512, 0.9900000000
+513, 0.9860000000
+514, 0.9810000000
+515, 0.9750000000
+516, 0.9680000000
+517, 0.9610000000
+518, 0.9530000000
+519, 0.9440000000
+520, 0.9350000000
+521, 0.9250000000
+522, 0.9150000000
+523, 0.9040000000
+524, 0.8920000000
+525, 0.8800000000
+526, 0.8670000000
+527, 0.8540000000
+528, 0.8400000000
+529, 0.8260000000
+530, 0.8110000000
+531, 0.7960000000
+532, 0.7810000000
+533, 0.7650000000
+534, 0.7490000000
+535, 0.7330000000
+536, 0.7170000000
+537, 0.7000000000
+538, 0.6830000000
+539, 0.6670000000
+540, 0.6500000000
+541, 0.6330000000
+542, 0.6160000000
+543, 0.5990000000
+544, 0.5810000000
+545, 0.5640000000
+546, 0.5480000000
+547, 0.5310000000
+548, 0.5140000000
+549, 0.4970000000
+550, 0.4810000000
+551, 0.4650000000
+552, 0.4480000000
+553, 0.4330000000
+554, 0.4170000000
+555, 0.4020000000
+556, 0.3864000000
+557, 0.3715000000
+558, 0.3569000000
+559, 0.3427000000
+560, 0.3288000000
+561, 0.3151000000
+562, 0.3018000000
+563, 0.2888000000
+564, 0.2762000000
+565, 0.2639000000
+566, 0.2519000000
+567, 0.2403000000
+568, 0.2291000000
+569, 0.2182000000
+570, 0.2076000000
+571, 0.1974000000
+572, 0.1876000000
+573, 0.1782000000
+574, 0.1690000000
+575, 0.1602000000
+576, 0.1517000000
+577, 0.1436000000
+578, 0.1358000000
+579, 0.1284000000
+580, 0.1212000000
+581, 0.1143000000
+582, 0.1078000000
+583, 0.1015000000
+584, 0.0956000000
+585, 0.0899000000
+586, 0.0845000000
+587, 0.0793000000
+588, 0.0745000000
+589, 0.0699000000
+590, 0.0655000000
+591, 0.0613000000
+592, 0.0574000000
+593, 0.0537000000
+594, 0.0502000000
+595, 0.0469000000
+596, 0.0438000000
+597, 0.0409000000
+598, 0.0381600000
+599, 0.0355800000
+600, 0.0331500000
+601, 0.0308700000
+602, 0.0287400000
+603, 0.0267400000
+604, 0.0248700000
+605, 0.0231200000
+606, 0.0214700000
+607, 0.0199400000
+608, 0.0185100000
+609, 0.0171800000
+610, 0.0159300000
+611, 0.0147700000
+612, 0.0136900000
+613, 0.0126900000
+614, 0.0117500000
+615, 0.0108800000
+616, 0.0100700000
+617, 0.0093200000
+618, 0.0086200000
+619, 0.0079700000
+620, 0.0073700000
+621, 0.0068200000
+622, 0.0063000000
+623, 0.0058200000
+624, 0.0053800000
+625, 0.0049700000
+626, 0.0045900000
+627, 0.0042400000
+628, 0.0039130000
+629, 0.0036130000
+630, 0.0033350000
+631, 0.0030790000
+632, 0.0028420000
+633, 0.0026230000
+634, 0.0024210000
+635, 0.0022350000
+636, 0.0020620000
+637, 0.0019030000
+638, 0.0017570000
+639, 0.0016210000
+640, 0.0014970000
+641, 0.0013820000
+642, 0.0012760000
+643, 0.0011780000
+644, 0.0010880000
+645, 0.0010050000
+646, 0.0009280000
+647, 0.0008570000
+648, 0.0007920000
+649, 0.0007320000
+650, 0.0006770000
+651, 0.0006260000
+652, 0.0005790000
+653, 0.0005360000
+654, 0.0004960000
+655, 0.0004590000
+656, 0.0004250000
+657, 0.0003935000
+658, 0.0003645000
+659, 0.0003377000
+660, 0.0003129000
+661, 0.0002901000
+662, 0.0002689000
+663, 0.0002493000
+664, 0.0002313000
+665, 0.0002146000
+666, 0.0001991000
+667, 0.0001848000
+668, 0.0001716000
+669, 0.0001593000
+670, 0.0001480000
+671, 0.0001375000
+672, 0.0001277000
+673, 0.0001187000
+674, 0.0001104000
+675, 0.0001026000
+676, 0.0000954000
+677, 0.0000888000
+678, 0.0000826000
+679, 0.0000769000
+680, 0.0000715000
+681, 0.0000666000
+682, 0.0000620000
+683, 0.0000578000
+684, 0.0000538000
+685, 0.0000501000
+686, 0.0000467000
+687, 0.0000436000
+688, 0.0000406000
+689, 0.0000378900
+690, 0.0000353300
+691, 0.0000329500
+692, 0.0000307500
+693, 0.0000287000
+694, 0.0000267900
+695, 0.0000250100
+696, 0.0000233600
+697, 0.0000218200
+698, 0.0000203800
+699, 0.0000190500
+700, 0.0000178000
+701, 0.0000166400
+702, 0.0000155600
+703, 0.0000145400
+704, 0.0000136000
+705, 0.0000127300
+706, 0.0000119100
+707, 0.0000111400
+708, 0.0000104300
+709, 0.0000097600
+710, 0.0000091400
+711, 0.0000085600
+712, 0.0000080200
+713, 0.0000075100
+714, 0.0000070400
+715, 0.0000066000
+716, 0.0000061800
+717, 0.0000058000
+718, 0.0000054400
+719, 0.0000051000
+720, 0.0000047800
+721, 0.0000044900
+722, 0.0000042100
+723, 0.0000039510
+724, 0.0000037090
+725, 0.0000034820
+726, 0.0000032700
+727, 0.0000030700
+728, 0.0000028840
+729, 0.0000027100
+730, 0.0000025460
+731, 0.0000023930
+732, 0.0000022500
+733, 0.0000021150
+734, 0.0000019890
+735, 0.0000018700
+736, 0.0000017590
+737, 0.0000016550
+738, 0.0000015570
+739, 0.0000014660
+740, 0.0000013790
+741, 0.0000012990
+742, 0.0000012230
+743, 0.0000011510
+744, 0.0000010840
+745, 0.0000010220
+746, 0.0000009620
+747, 0.0000009070
+748, 0.0000008550
+749, 0.0000008060
+750, 0.0000007600
+751, 0.0000007160
+752, 0.0000006750
+753, 0.0000006370
+754, 0.0000006010
+755, 0.0000005670
+756, 0.0000005350
+757, 0.0000005050
+758, 0.0000004770
+759, 0.0000004500
+760, 0.0000004250
+761, 0.0000004010
+762, 0.0000003790
+763, 0.0000003580
+764, 0.0000003382
+765, 0.0000003196
+766, 0.0000003021
+767, 0.0000002855
+768, 0.0000002699
+769, 0.0000002552
+770, 0.0000002413
+771, 0.0000002282
+772, 0.0000002159
+773, 0.0000002042
+774, 0.0000001932
+775, 0.0000001829
+776, 0.0000001731
+777, 0.0000001638
+778, 0.0000001551
+779, 0.0000001468
+780, 0.0000001390 \ No newline at end of file
diff --git a/color/vl1924e_1.csv b/color/vl1924e_1.csv
new file mode 100644
index 0000000..8b240c0
--- /dev/null
+++ b/color/vl1924e_1.csv
@@ -0,0 +1,471 @@
+360, 0.0000039170000
+361, 0.0000043935810
+362, 0.0000049296040
+363, 0.0000055321360
+364, 0.0000062082450
+365, 0.0000069650000
+366, 0.0000078132190
+367, 0.0000087673360
+368, 0.0000098398440
+369, 0.0000110432300
+370, 0.0000123900000
+371, 0.0000138864100
+372, 0.0000155572800
+373, 0.0000174429600
+374, 0.0000195837500
+375, 0.0000220200000
+376, 0.0000248396500
+377, 0.0000280412600
+378, 0.0000315310400
+379, 0.0000352152100
+380, 0.0000390000000
+381, 0.0000428264000
+382, 0.0000469146000
+383, 0.0000515896000
+384, 0.0000571764000
+385, 0.0000640000000
+386, 0.0000723442100
+387, 0.0000822122400
+388, 0.0000935081600
+389, 0.0001061361000
+390, 0.0001200000000
+391, 0.0001349840000
+392, 0.0001514920000
+393, 0.0001702080000
+394, 0.0001918160000
+395, 0.0002170000000
+396, 0.0002469067000
+397, 0.0002812400000
+398, 0.0003185200000
+399, 0.0003572667000
+400, 0.0003960000000
+401, 0.0004337147000
+402, 0.0004730240000
+403, 0.0005178760000
+404, 0.0005722187000
+405, 0.0006400000000
+406, 0.0007245600000
+407, 0.0008255000000
+408, 0.0009411600000
+409, 0.0010698800000
+410, 0.0012100000000
+411, 0.0013620910000
+412, 0.0015307520000
+413, 0.0017203680000
+414, 0.0019353230000
+415, 0.0021800000000
+416, 0.0024548000000
+417, 0.0027640000000
+418, 0.0031178000000
+419, 0.0035264000000
+420, 0.0040000000000
+421, 0.0045462400000
+422, 0.0051593200000
+423, 0.0058292800000
+424, 0.0065461600000
+425, 0.0073000000000
+426, 0.0080865070000
+427, 0.0089087200000
+428, 0.0097676800000
+429, 0.0106644300000
+430, 0.0116000000000
+431, 0.0125731700000
+432, 0.0135827200000
+433, 0.0146296800000
+434, 0.0157150900000
+435, 0.0168400000000
+436, 0.0180073600000
+437, 0.0192144800000
+438, 0.0204539200000
+439, 0.0217182400000
+440, 0.0230000000000
+441, 0.0242946100000
+442, 0.0256102400000
+443, 0.0269585700000
+444, 0.0283512500000
+445, 0.0298000000000
+446, 0.0313108300000
+447, 0.0328836800000
+448, 0.0345211200000
+449, 0.0362257100000
+450, 0.0380000000000
+451, 0.0398466700000
+452, 0.0417680000000
+453, 0.0437660000000
+454, 0.0458426700000
+455, 0.0480000000000
+456, 0.0502436800000
+457, 0.0525730400000
+458, 0.0549805600000
+459, 0.0574587200000
+460, 0.0600000000000
+461, 0.0626019700000
+462, 0.0652775200000
+463, 0.0680420800000
+464, 0.0709110900000
+465, 0.0739000000000
+466, 0.0770160000000
+467, 0.0802664000000
+468, 0.0836668000000
+469, 0.0872328000000
+470, 0.0909800000000
+471, 0.0949175500000
+472, 0.0990458400000
+473, 0.1033674000000
+474, 0.1078846000000
+475, 0.1126000000000
+476, 0.1175320000000
+477, 0.1226744000000
+478, 0.1279928000000
+479, 0.1334528000000
+480, 0.1390200000000
+481, 0.1446764000000
+482, 0.1504693000000
+483, 0.1564619000000
+484, 0.1627177000000
+485, 0.1693000000000
+486, 0.1762431000000
+487, 0.1835581000000
+488, 0.1912735000000
+489, 0.1994180000000
+490, 0.2080200000000
+491, 0.2171199000000
+492, 0.2267345000000
+493, 0.2368571000000
+494, 0.2474812000000
+495, 0.2586000000000
+496, 0.2701849000000
+497, 0.2822939000000
+498, 0.2950505000000
+499, 0.3085780000000
+500, 0.3230000000000
+501, 0.3384021000000
+502, 0.3546858000000
+503, 0.3716986000000
+504, 0.3892875000000
+505, 0.4073000000000
+506, 0.4256299000000
+507, 0.4443096000000
+508, 0.4633944000000
+509, 0.4829395000000
+510, 0.5030000000000
+511, 0.5235693000000
+512, 0.5445120000000
+513, 0.5656900000000
+514, 0.5869653000000
+515, 0.6082000000000
+516, 0.6293456000000
+517, 0.6503068000000
+518, 0.6708752000000
+519, 0.6908424000000
+520, 0.7100000000000
+521, 0.7281852000000
+522, 0.7454636000000
+523, 0.7619694000000
+524, 0.7778368000000
+525, 0.7932000000000
+526, 0.8081104000000
+527, 0.8224962000000
+528, 0.8363068000000
+529, 0.8494916000000
+530, 0.8620000000000
+531, 0.8738108000000
+532, 0.8849624000000
+533, 0.8954936000000
+534, 0.9054432000000
+535, 0.9148501000000
+536, 0.9237348000000
+537, 0.9320924000000
+538, 0.9399226000000
+539, 0.9472252000000
+540, 0.9540000000000
+541, 0.9602561000000
+542, 0.9660074000000
+543, 0.9712606000000
+544, 0.9760225000000
+545, 0.9803000000000
+546, 0.9840924000000
+547, 0.9874182000000
+548, 0.9903128000000
+549, 0.9928116000000
+550, 0.9949501000000
+551, 0.9967108000000
+552, 0.9980983000000
+553, 0.9991120000000
+554, 0.9997482000000
+555, 1.0000000000000
+556, 0.9998567000000
+557, 0.9993046000000
+558, 0.9983255000000
+559, 0.9968987000000
+560, 0.9950000000000
+561, 0.9926005000000
+562, 0.9897426000000
+563, 0.9864444000000
+564, 0.9827241000000
+565, 0.9786000000000
+566, 0.9740837000000
+567, 0.9691712000000
+568, 0.9638568000000
+569, 0.9581349000000
+570, 0.9520000000000
+571, 0.9454504000000
+572, 0.9384992000000
+573, 0.9311628000000
+574, 0.9234576000000
+575, 0.9154000000000
+576, 0.9070064000000
+577, 0.8982772000000
+578, 0.8892048000000
+579, 0.8797816000000
+580, 0.8700000000000
+581, 0.8598613000000
+582, 0.8493920000000
+583, 0.8386220000000
+584, 0.8275813000000
+585, 0.8163000000000
+586, 0.8047947000000
+587, 0.7930820000000
+588, 0.7811920000000
+589, 0.7691547000000
+590, 0.7570000000000
+591, 0.7447541000000
+592, 0.7324224000000
+593, 0.7200036000000
+594, 0.7074965000000
+595, 0.6949000000000
+596, 0.6822192000000
+597, 0.6694716000000
+598, 0.6566744000000
+599, 0.6438448000000
+600, 0.6310000000000
+601, 0.6181555000000
+602, 0.6053144000000
+603, 0.5924756000000
+604, 0.5796379000000
+605, 0.5668000000000
+606, 0.5539611000000
+607, 0.5411372000000
+608, 0.5283528000000
+609, 0.5156323000000
+610, 0.5030000000000
+611, 0.4904688000000
+612, 0.4780304000000
+613, 0.4656776000000
+614, 0.4534032000000
+615, 0.4412000000000
+616, 0.4290800000000
+617, 0.4170360000000
+618, 0.4050320000000
+619, 0.3930320000000
+620, 0.3810000000000
+621, 0.3689184000000
+622, 0.3568272000000
+623, 0.3447768000000
+624, 0.3328176000000
+625, 0.3210000000000
+626, 0.3093381000000
+627, 0.2978504000000
+628, 0.2865936000000
+629, 0.2756245000000
+630, 0.2650000000000
+631, 0.2547632000000
+632, 0.2448896000000
+633, 0.2353344000000
+634, 0.2260528000000
+635, 0.2170000000000
+636, 0.2081616000000
+637, 0.1995488000000
+638, 0.1911552000000
+639, 0.1829744000000
+640, 0.1750000000000
+641, 0.1672235000000
+642, 0.1596464000000
+643, 0.1522776000000
+644, 0.1451259000000
+645, 0.1382000000000
+646, 0.1315003000000
+647, 0.1250248000000
+648, 0.1187792000000
+649, 0.1127691000000
+650, 0.1070000000000
+651, 0.1014762000000
+652, 0.0961886400000
+653, 0.0911229600000
+654, 0.0862648500000
+655, 0.0816000000000
+656, 0.0771206400000
+657, 0.0728255200000
+658, 0.0687100800000
+659, 0.0647697600000
+660, 0.0610000000000
+661, 0.0573962100000
+662, 0.0539550400000
+663, 0.0506737600000
+664, 0.0475496500000
+665, 0.0445800000000
+666, 0.0417587200000
+667, 0.0390849600000
+668, 0.0365638400000
+669, 0.0342004800000
+670, 0.0320000000000
+671, 0.0299626100000
+672, 0.0280766400000
+673, 0.0263293600000
+674, 0.0247080500000
+675, 0.0232000000000
+676, 0.0218007700000
+677, 0.0205011200000
+678, 0.0192810800000
+679, 0.0181206900000
+680, 0.0170000000000
+681, 0.0159037900000
+682, 0.0148371800000
+683, 0.0138106800000
+684, 0.0128347800000
+685, 0.0119200000000
+686, 0.0110683100000
+687, 0.0102733900000
+688, 0.0095333110000
+689, 0.0088461570000
+690, 0.0082100000000
+691, 0.0076237810000
+692, 0.0070854240000
+693, 0.0065914760000
+694, 0.0061384850000
+695, 0.0057230000000
+696, 0.0053430590000
+697, 0.0049957960000
+698, 0.0046764040000
+699, 0.0043800750000
+700, 0.0041020000000
+701, 0.0038384530000
+702, 0.0035890990000
+703, 0.0033542190000
+704, 0.0031340930000
+705, 0.0029290000000
+706, 0.0027381390000
+707, 0.0025598760000
+708, 0.0023932440000
+709, 0.0022372750000
+710, 0.0020910000000
+711, 0.0019535870000
+712, 0.0018245800000
+713, 0.0017035800000
+714, 0.0015901870000
+715, 0.0014840000000
+716, 0.0013844960000
+717, 0.0012912680000
+718, 0.0012040920000
+719, 0.0011227440000
+720, 0.0010470000000
+721, 0.0009765896000
+722, 0.0009111088000
+723, 0.0008501332000
+724, 0.0007932384000
+725, 0.0007400000000
+726, 0.0006900827000
+727, 0.0006433100000
+728, 0.0005994960000
+729, 0.0005584547000
+730, 0.0005200000000
+731, 0.0004839136000
+732, 0.0004500528000
+733, 0.0004183452000
+734, 0.0003887184000
+735, 0.0003611000000
+736, 0.0003353835000
+737, 0.0003114404000
+738, 0.0002891656000
+739, 0.0002684539000
+740, 0.0002492000000
+741, 0.0002313019000
+742, 0.0002146856000
+743, 0.0001992884000
+744, 0.0001850475000
+745, 0.0001719000000
+746, 0.0001597781000
+747, 0.0001486044000
+748, 0.0001383016000
+749, 0.0001287925000
+750, 0.0001200000000
+751, 0.0001118595000
+752, 0.0001043224000
+753, 0.0000973356000
+754, 0.0000908458700
+755, 0.0000848000000
+756, 0.0000791466700
+757, 0.0000738580000
+758, 0.0000689160000
+759, 0.0000643026700
+760, 0.0000600000000
+761, 0.0000559818700
+762, 0.0000522256000
+763, 0.0000487184000
+764, 0.0000454474700
+765, 0.0000424000000
+766, 0.0000395610400
+767, 0.0000369151200
+768, 0.0000344486800
+769, 0.0000321481600
+770, 0.0000300000000
+771, 0.0000279912500
+772, 0.0000261135600
+773, 0.0000243602400
+774, 0.0000227246100
+775, 0.0000212000000
+776, 0.0000197785500
+777, 0.0000184528500
+778, 0.0000172168700
+779, 0.0000160645900
+780, 0.0000149900000
+781, 0.0000139872800
+782, 0.0000130515500
+783, 0.0000121781800
+784, 0.0000113625400
+785, 0.0000106000000
+786, 0.0000098858770
+787, 0.0000092173040
+788, 0.0000085923620
+789, 0.0000080091330
+790, 0.0000074657000
+791, 0.0000069595670
+792, 0.0000064879950
+793, 0.0000060486990
+794, 0.0000056393960
+795, 0.0000052578000
+796, 0.0000049017710
+797, 0.0000045697200
+798, 0.0000042601940
+799, 0.0000039717390
+800, 0.0000037029000
+801, 0.0000034521630
+802, 0.0000032183020
+803, 0.0000030003000
+804, 0.0000027971390
+805, 0.0000026078000
+806, 0.0000024312200
+807, 0.0000022665310
+808, 0.0000021130130
+809, 0.0000019699430
+810, 0.0000018366000
+811, 0.0000017122300
+812, 0.0000015962280
+813, 0.0000014880900
+814, 0.0000013873140
+815, 0.0000012934000
+816, 0.0000012058200
+817, 0.0000011241430
+818, 0.0000010480090
+819, 0.0000009770578
+820, 0.0000009109300
+821, 0.0000008492513
+822, 0.0000007917212
+823, 0.0000007380904
+824, 0.0000006881098
+825, 0.0000006415300
+826, 0.0000005980895
+827, 0.0000005575746
+828, 0.0000005198080
+829, 0.0000004846123
+830, 0.0000004518100 \ No newline at end of file
diff --git a/detectors/__init__.py b/detectors/__init__.py
index 147fd59..d7583a6 100644
--- a/detectors/__init__.py
+++ b/detectors/__init__.py
@@ -24,6 +24,10 @@ def build_lbne_200kton():
def build_minilbne():
return build_lbne(radius/10, height/10, nstrings//10, pmts_per_string//10, endcap_spacing)
+@buildable('microlbne')
+def build_microlbne():
+ return build_lbne(radius/40, height/40, nstrings//40, pmts_per_string//40, endcap_spacing)
+
@buildable('sno')
def build_sno():
return build_sno_detector()
diff --git a/geometry.py b/geometry.py
index 1709a87..5fdfa1f 100644
--- a/geometry.py
+++ b/geometry.py
@@ -4,6 +4,7 @@ import pycuda.driver as cuda
from pycuda import gpuarray
from itertools import *
from itertoolset import *
+from camera import timeit
# all material/surface properties are interpolated at these
# wavelengths when they are sent to the gpu
@@ -213,6 +214,7 @@ class Geometry(object):
return len(self.solids)-1
+ @timeit
def build(self, bits=8, shift=3):
offsets = [ (0,0) ]
for solid in self.solids:
@@ -352,18 +354,6 @@ class Geometry(object):
if unique_bit_shifted_zvalues.size == 1:
break
- def load_colors(self, module):
- self.colors = \
- np.concatenate([solid.color for solid in self.solids])[reorder]
-
- if not hasattr(self, 'colors_gpu'):
- raise Exception('cannot load colors before geometry is loaded.')
-
- self.colors_gpu.set(self.colors.astype(np.uint32))
-
- set_colors = module.get_function('set_colors')
- set_colors(self.colors_gpu, block=(1,1,1), grid=(1,1))
-
def load(self, module):
"""
Load the bounding volume hierarchy onto the GPU module,
diff --git a/solids/__init__.py b/solids/__init__.py
index 9de5c8d..479bacd 100644
--- a/solids/__init__.py
+++ b/solids/__init__.py
@@ -27,6 +27,12 @@ def build_12inch_pmt_with_lc(outer_material=water, theta=np.pi/8):
pmt = build_12inch_pmt(outer_material, theta)
return pmt + build_light_collector(pmt, a=lc_12inch_a, b=lc_12inch_b, d=lc_12inch_d, rmin=lc_12inch_rmin, rmax=lc_12inch_rmax)
+@buildable('12inch_pmt_with_lc_hd')
+def build_12inch_pmt_with_lc_hd(outer_material=water, theta=np.pi/64):
+ pmt = build_12inch_pmt(outer_material, theta)
+ return pmt + build_light_collector(pmt, a=lc_12inch_a, b=lc_12inch_b, d=lc_12inch_d, rmin=lc_12inch_rmin, rmax=lc_12inch_rmax, npoints=100)
+
+
@buildable('8inch_pmt')
def build_8inch_pmt(outer_material=water, theta=np.pi/12):
return build_pmt(dir + '/sno_pmt_reduced.txt', 0.003, outer_material, theta)
diff --git a/src/kernel.cu b/src/kernel.cu
index a36ff91..82efe88 100644
--- a/src/kernel.cu
+++ b/src/kernel.cu
@@ -9,14 +9,12 @@
#include "photon.h"
#include "alpha.h"
-#define RED_WAVELENGTH 685
-#define BLUE_WAVELENGTH 445
-#define GREEN_WAVELENGTH 545
-
__device__ void fAtomicAdd(float *addr, float data)
{
while (data)
+ {
data = atomicExch(addr, data+atomicExch(addr, 0.0f));
+ }
}
__device__ __noinline__ void to_diffuse(Photon &p, State &s, curandState &rng, const int &max_steps)
@@ -28,9 +26,9 @@ __device__ __noinline__ void to_diffuse(Photon &p, State &s, curandState &rng, c
int command;
- command = fill_state(s, p, rng);
+ fill_state(s, p);
- if (command == BREAK)
+ if (p.last_hit_triangle == -1)
break;
command = propagate_to_boundary(p, s, rng);
@@ -64,127 +62,111 @@ __device__ __noinline__ void to_diffuse(Photon &p, State &s, curandState &rng, c
extern "C"
{
-/* Translate `points` by the vector `v` */
-__global__ void translate(int nthreads, float3 *points, float3 v)
+/* Translate the points `a` by the vector `v` */
+__global__ void translate(int nthreads, float3 *a, float3 v)
{
int id = blockIdx.x*blockDim.x + threadIdx.x;
if (id >= nthreads)
return;
- points[id] += v;
+ a[id] += v;
}
-/* Rotate `points` through an angle `phi` counter-clockwise about the
+/* Rotate the points `a` through an angle `phi` counter-clockwise about the
axis `axis` (when looking towards +infinity). */
-__global__ void rotate(int nthreads, float3 *points, float phi, float3 axis)
+__global__ void rotate(int nthreads, float3 *a, float phi, float3 axis)
{
int id = blockIdx.x*blockDim.x + threadIdx.x;
if (id >= nthreads)
return;
- points[id] = rotate(points[id], phi, axis);
+ a[id] = rotate(a[id], phi, axis);
}
- __global__ void build_rgb_lookup(int nthreads, float3 *positions, float3 *directions, float3 *rgb_lookup1, float3 *rgb_lookup2, unsigned long long seed, int runs, int max_steps)
+__global__ void rotate_around_point(int nthreads, float3 *a, float phi, float3 axis, float3 point)
{
int id = blockIdx.x*blockDim.x + threadIdx.x;
if (id >= nthreads)
return;
- curandState rng;
- curand_init(seed, id, 0, &rng);
-
- Photon p0;
- p0.position = positions[id];
- p0.direction = directions[id];
- p0.direction /= norm(p0.direction);
- p0.polarization = uniform_sphere(&rng);
- p0.time = 0.0f;
- p0.history = 0x0;
-
- float distance;
+ a[id] -= point;
+ a[id] = rotate(a[id], phi, axis);
+ a[id] += point;
+}
- int hit_triangle = intersect_mesh(p0.position, p0.direction, distance);
+__global__ void update_xyz_lookup(int nthreads, int total_threads, int offset, float3 position, curandState *rng_states, float wavelength, float3 xyz, float3 *xyz_lookup1, float3 *xyz_lookup2, int max_steps)
+{
+ int kernel_id = blockIdx.x*blockDim.x + threadIdx.x;
+ int id = kernel_id + offset;
- if (hit_triangle != id)
+ if (kernel_id >= nthreads || id >= total_threads)
return;
- // note triangles from built geometry mesh are in order
+ curandState rng = rng_states[kernel_id];
- uint4 triangle_data = g_triangles[hit_triangle];
+ uint4 triangle_data = g_triangles[id];
float3 v0 = g_vertices[triangle_data.x];
float3 v1 = g_vertices[triangle_data.y];
float3 v2 = g_vertices[triangle_data.z];
- float cos_theta = dot(normalize(cross(v1-v0, v2-v1)), -p0.direction);
+ float a = curand_uniform(&rng);
+ float b = uniform(&rng, 0.0f, (1.0f - a));
+ float c = 1.0f - a - b;
- if (cos_theta < 0.0f)
- cos_theta = dot(-normalize(cross(v1-v0, v2-v1)), -p0.direction);
+ float3 direction = a*v0 + b*v1 + c*v2 - position;
+ direction /= norm(direction);
- Photon p;
- State s;
+ float distance;
+ int triangle_index = intersect_mesh(position, direction, distance);
- for (int i=0; i < runs; i++)
+ if (triangle_index != id)
{
- p = p0;
- p.wavelength = RED_WAVELENGTH;
+ rng_states[kernel_id] = rng;
+ return;
+ }
- to_diffuse(p, s, rng, max_steps);
+ float cos_theta = dot(normalize(cross(v1-v0,v2-v1)),-direction);
- if (p.history & REFLECT_DIFFUSE)
- {
- if (s.inside_to_outside)
- {
- fAtomicAdd(&rgb_lookup1[p.last_hit_triangle].x, cos_theta);
- }
- else
- {
- fAtomicAdd(&rgb_lookup2[p.last_hit_triangle].x, cos_theta);
- }
- }
+ if (cos_theta < 0.0f)
+ cos_theta = dot(-normalize(cross(v1-v0,v2-v1)),-direction);
- p = p0;
- p.wavelength = GREEN_WAVELENGTH;
+ Photon p;
+ p.position = position;
+ p.direction = direction;
+ p.wavelength = wavelength;
+ p.polarization = uniform_sphere(&rng);
+ p.last_hit_triangle = -1;
+ p.time = 0;
+ p.history = 0;
- to_diffuse(p, s, rng, max_steps);
+ State s;
+ to_diffuse(p, s, rng, max_steps);
- if (p.history & REFLECT_DIFFUSE)
+ if (p.history & REFLECT_DIFFUSE)
+ {
+ if (s.inside_to_outside)
{
- if (s.inside_to_outside)
- {
- fAtomicAdd(&rgb_lookup1[p.last_hit_triangle].y, cos_theta);
- }
- else
- {
- fAtomicAdd(&rgb_lookup2[p.last_hit_triangle].y, cos_theta);
- }
+ fAtomicAdd(&xyz_lookup1[p.last_hit_triangle].x, cos_theta*xyz.x);
+ fAtomicAdd(&xyz_lookup1[p.last_hit_triangle].y, cos_theta*xyz.y);
+ fAtomicAdd(&xyz_lookup1[p.last_hit_triangle].z, cos_theta*xyz.z);
}
-
- p = p0;
- p.wavelength = BLUE_WAVELENGTH;
-
- to_diffuse(p, s, rng, max_steps);
-
- if (p.history & REFLECT_DIFFUSE)
+ else
{
- if (s.inside_to_outside)
- {
- fAtomicAdd(&rgb_lookup1[p.last_hit_triangle].z, cos_theta);
- }
- else
- {
- fAtomicAdd(&rgb_lookup2[p.last_hit_triangle].z, cos_theta);
- }
+ fAtomicAdd(&xyz_lookup2[p.last_hit_triangle].x, cos_theta*xyz.x);
+ fAtomicAdd(&xyz_lookup2[p.last_hit_triangle].y, cos_theta*xyz.y);
+ fAtomicAdd(&xyz_lookup2[p.last_hit_triangle].z, cos_theta*xyz.z);
}
}
-} // build_rgb_lookup
+ rng_states[kernel_id] = rng;
+
+} // update_xyz_lookup
-__global__ void render(int nthreads, curandState *rng_states, float3 *positions, float3 *directions, float *red, float *green, float *blue, float3 *rgb_lookup1, float3 *rgb_lookup2, int max_steps)
+__global__ void update_xyz_image(int nthreads, curandState *rng_states, float3 *positions, float3 *directions, float wavelength, float3 xyz, float3 *xyz_lookup1, float3 *xyz_lookup2, float3 *image, int nlookup_calls, int max_steps)
{
int id = blockIdx.x*blockDim.x + threadIdx.x;
@@ -193,81 +175,69 @@ __global__ void render(int nthreads, curandState *rng_states, float3 *positions,
curandState rng = rng_states[id];
- Photon p0;
- p0.position = positions[id];
- p0.direction = directions[id];
- p0.direction /= norm(p0.direction);
- p0.polarization = uniform_sphere(&rng);
- p0.time = 0.0f;
- p0.history = 0x0;
+ Photon p;
+ p.position = positions[id];
+ p.direction = directions[id];
+ p.direction /= norm(p.direction);
+ p.wavelength = wavelength;
+ p.polarization = uniform_sphere(&rng);
+ p.last_hit_triangle = -1;
+ p.time = 0;
+ p.history = 0;
State s;
- Photon p = p0;
- p.wavelength = RED_WAVELENGTH;
-
to_diffuse(p, s, rng, max_steps);
if (p.history & REFLECT_DIFFUSE)
{
if (s.inside_to_outside)
{
- red[id] = rgb_lookup1[p.last_hit_triangle].x;
- }
- else
- {
- red[id] = rgb_lookup2[p.last_hit_triangle].x;
- }
- }
- else
- {
- red[id] = 0.0;
- }
-
- p = p0;
- p.wavelength = GREEN_WAVELENGTH;
-
- to_diffuse(p, s, rng, max_steps);
-
- if (p.history & REFLECT_DIFFUSE)
- {
- if (s.inside_to_outside)
- {
- green[id] = rgb_lookup1[p.last_hit_triangle].y;
+ image[id].x += xyz.x*xyz_lookup1[p.last_hit_triangle].x/nlookup_calls;
+ image[id].y += xyz.y*xyz_lookup1[p.last_hit_triangle].y/nlookup_calls;
+ image[id].z += xyz.z*xyz_lookup1[p.last_hit_triangle].z/nlookup_calls;
}
else
{
- green[id] = rgb_lookup2[p.last_hit_triangle].y;
+ image[id].x += xyz.x*xyz_lookup2[p.last_hit_triangle].x/nlookup_calls;
+ image[id].y += xyz.y*xyz_lookup2[p.last_hit_triangle].y/nlookup_calls;
+ image[id].z += xyz.z*xyz_lookup2[p.last_hit_triangle].z/nlookup_calls;
}
}
- else
- {
- green[id] = 0.0;
- }
-
- p = p0;
- p.wavelength = BLUE_WAVELENGTH;
-
- to_diffuse(p, s, rng, max_steps);
-
- if (p.history & REFLECT_DIFFUSE)
- {
- if (s.inside_to_outside)
- {
- blue[id] = rgb_lookup1[p.last_hit_triangle].z;
- }
- else
- {
- blue[id] = rgb_lookup2[p.last_hit_triangle].z;
- }
- }
- else
- {
- blue[id] = 0.0;
- }
-
+
rng_states[id] = rng;
-} // render
+} // update_xyz_image
+
+__global__ void process_image(int nthreads, float3 *image, int *pixels, int nimages)
+{
+ int id = blockIdx.x*blockDim.x + threadIdx.x;
+
+ if (id >= nthreads)
+ return;
+
+ float3 rgb = image[id]/nimages;
+
+ if (rgb.x < 0.0f)
+ rgb.x = 0.0f;
+ if (rgb.y < 0.0f)
+ rgb.y = 0.0f;
+ if (rgb.z < 0.0f)
+ rgb.z = 0.0f;
+
+ if (rgb.x > 1.0f)
+ rgb.x = 1.0f;
+ if (rgb.y > 1.0f)
+ rgb.y = 1.0f;
+ if (rgb.z > 1.0f)
+ rgb.z = 1.0f;
+
+ unsigned int r = floorf(rgb.x*255.0f);
+ unsigned int g = floorf(rgb.y*255.0f);
+ unsigned int b = floorf(rgb.z*255.0f);
+
+ pixels[id] = r << 16 | g << 8 | b;
+
+} // process_image
/* Trace the rays starting at `positions` traveling in the direction
`directions` to their intersection with the global mesh. If the ray
@@ -322,9 +292,9 @@ __global__ void propagate(int nthreads, curandState *rng_states, float3 *positio
int command;
- command = fill_state(s, p, rng);
+ fill_state(s, p);
- if (command == BREAK)
+ if (p.last_hit_triangle == -1)
break;
command = propagate_to_boundary(p, s, rng);
diff --git a/src/photon.h b/src/photon.h
index ad4c26c..f471866 100644
--- a/src/photon.h
+++ b/src/photon.h
@@ -19,8 +19,6 @@ struct Photon
unsigned int history;
int last_hit_triangle;
-
- //curandState rng;
};
struct State
@@ -51,14 +49,14 @@ enum
enum {BREAK, CONTINUE, PASS}; // return value from propagate_to_boundary
-__device__ int fill_state(State &s, Photon &p, curandState &rng)
+__device__ void fill_state(State &s, Photon &p)
{
p.last_hit_triangle = intersect_mesh(p.position, p.direction, s.distance_to_boundary, p.last_hit_triangle);
if (p.last_hit_triangle == -1)
{
p.history |= NO_HIT;
- return BREAK;
+ return;
}
uint4 triangle_data = g_triangles[p.last_hit_triangle];
@@ -98,8 +96,6 @@ __device__ int fill_state(State &s, Photon &p, curandState &rng)
s.absorption_length = interp_property(p.wavelength, material1.absorption_length);
s.scattering_length = interp_property(p.wavelength, material1.scattering_length);
- return PASS;
-
} // fill_state
__device__ void rayleigh_scatter(Photon &p, curandState &rng)
generated by cgit (git 2.34.1) at 2025-05-08 07:56:37 +0000