added a camera class which is able to render using the simple ray tracer or the hybrid monte carlo ray tracer in a separate thread. the camera object is initialized by passing a pycuda.driver.Context object and a threading.Lock object; you can then run kernels and copy arrays to and from the device in the same context as that used by the camera by acquiring the lock, calling context.push(), executing the kernel and/or copying arrays, calling context.pop(), and then releasing the lock. fixed mistakes in both build_rgb_lookup() and render() where I had accidently switched the green and blue photons. updated the lbne geometry with the latest specifications. added profiles for the sno pmt and concentrator.

10 files changed, 551 insertions, 42 deletions

diff --git a/camera.py b/camera.py
index 5a05535..57d0a67 100644
--- a/camera.py
+++ b/camera.py
@@ -1,11 +1,25 @@
 import numpy as np
-from itertools import product
+from itertools import product, count
+from threading import Thread, Lock
+import time
+import os
 
-class Camera(object):
+from transform import rotate
+
+import pygame
+from pygame.locals import *
+
+from pycuda import gpuarray
+from pycuda.characterize import sizeof
+import pycuda.driver as cuda
+
+def get_rays(position, size = (800, 600), film_size = (0.035, 0.024), focal_length=0.05):
     """
-    Pinhole camera object.
+    Generate ray positions and directions from a pinhole camera facing the negative y direction.
 
     Args:
+        - position: tuple,
+            Position of the camera.
         - size: tuple, *optional*
             Pixel array shape.
         - film_size: tuple, *optional*
@@ -13,21 +27,363 @@ class Camera(object):
         - focal_length: float, *optional*
             Focal length of camera.
     """
-    def __init__(self, size = (800, 600), film_size = (0.035, 0.024), \
-                     focal_length=0.05):
-        x = np.linspace(-film_size[0]/2, film_size[0]/2, size[0])
-        z = np.linspace(-film_size[1]/2, film_size[1]/2, size[1])
 
-        self.grid = np.array(tuple(product(x,[0],z)))
+    x = np.linspace(-film_size[0]/2, film_size[0]/2, size[0])
+    z = np.linspace(-film_size[1]/2, film_size[1]/2, size[1])
+
+    grid = np.array(tuple(product(x,[0],z)))
+
+    grid += (0,focal_length,0)
+    focal_point = np.zeros(3)
+
+    grid += position
+    focal_point += position
+
+    return grid, focal_point-grid
+
+class Camera(Thread):
+    def __init__(self, geometry, module, context, lock, size=(800,600)):
+        Thread.__init__(self)
+        self.geometry = geometry
+        self.module = module
+        self.context = context
+        self.lock = lock
+
+        self.size = size
+        self.width, self.height = size
+
+        pygame.init()
+        self.screen = pygame.display.set_mode(size)
+        pygame.display.set_caption('')
+
+        with self.lock:
+            self.context.push()
+            self.ray_trace_kernel = self.module.get_function('ray_trace')
+            self.rotate_kernel = self.module.get_function('rotate')
+            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.init_rng_kernel = self.module.get_function('init_rng')
+            self.context.pop()
+
+        lower_bound, upper_bound = self.geometry.mesh.get_bounds()
+        self.scale = np.linalg.norm(upper_bound-lower_bound)
+        self.source_position = [0,0,upper_bound[2]+1.0]
+
+        self.nblocks = 64
+
+        with self.lock:
+            self.context.push()
+            print 'initializing random number states.'
+            rng_state_count = max(len(geometry.mesh.triangles), self.width*self.height)
+            self.rng_states_gpu = cuda.mem_alloc(rng_state_count*sizeof('curandStateXORWOW', '#include <curand_kernel.h>'))
+            self.init_rng_kernel(np.int32(rng_state_count), self.rng_states_gpu, np.int32(0), np.int32(0), block=(self.nblocks,1,1), grid=(rng_state_count//self.nblocks+1,1))
+            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(*self.source_position))
+
+            source_directions = np.mean(self.geometry.mesh[:], axis=1) - self.source_position
+            self.source_directions_gpu = gpuarray.to_gpu(source_directions.astype(np.float32).view(gpuarray.vec.float3))
+
+            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.max_steps = 10
+            rgb_runs = 100
+            
+            print 'building rgb lookup.'
+            self.build_rgb_lookup_kernel(np.int32(self.source_positions_gpu.size), self.rng_states_gpu, self.source_positions_gpu, self.source_directions_gpu, self.rgb_lookup1_gpu, self.rgb_lookup2_gpu, 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.'
+
+            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))
+
+            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))
+
+            self.nruns = 1
+
+            self.context.pop()
+
+        self.point = np.array([0, self.scale*1.75, (lower_bound[2]+upper_bound[2])/2])
+        self.axis1 = np.array([1,0,0], float)
+        self.axis2 = np.array([0,0,1], float)
+
+        origins, directions = get_rays(self.point)
+
+        with self.lock:
+            self.context.push()
+            self.origins_gpu = gpuarray.to_gpu(origins.astype(np.float32).view(gpuarray.vec.float3))
+            self.directions_gpu = gpuarray.to_gpu(directions.astype(np.float32).view(gpuarray.vec.float3))
+            self.pixels_gpu = gpuarray.zeros(self.width*self.height, dtype=np.int32)
+            self.context.pop()
+
+        self.render = False
+
+        self.render_run_count = 1
+
+    def screenshot(self, dir=''):
+        root, ext = 'screenshot', 'png'
+
+        for i in count():
+            filename = os.path.join(dir, '.'.join([root + str(i).zfill(4), ext]))
+
+            if not os.path.exists(filename):
+                break
+
+        pygame.image.save(self.screen, filename)
+        print 'image saved to %s' % filename
+
+    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.point = rotate(self.point, phi, n)
+            self.axis1 = rotate(self.axis1, phi, n)
+            self.axis2 = rotate(self.axis2, phi, n)
+            self.context.pop()
+
+    def translate(self, v):
+        with self.lock:
+            self.context.push()
+            self.translate_kernel(np.int32(self.pixels_gpu.size), self.origins_gpu, gpuarray.vec.make_float3(*v), block=(self.nblocks,1,1), grid=(self.pixels_gpu.size//self.nblocks,1))
+
+            self.point += v
+            self.context.pop()
+
+    def update(self):
+        with self.lock:
+            self.context.push()
+            t0 = time.time()
+            if self.render:
+                self.render_kernel(np.int32(self.pixels_gpu.size), self.rng_states_gpu, self.origins_gpu, self.directions_gpu, self.rgb_lookup1_gpu, self.rgb_lookup2_gpu, np.int32(self.render_run_count), self.pixels_gpu, np.int32(self.max_steps), block=(self.nblocks,1,1), grid=(self.pixels_gpu.size//self.nblocks+1,1))
+            else:
+                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 'elapsed %f sec.' % elapsed
+
+            pygame.surfarray.blit_array(self.screen, self.pixels_gpu.get().reshape(self.size))
+            pygame.display.flip()
+            self.context.pop()
+
+    def run(self):
+        self.update()
+        
+        done = False
+        clicked = False
+        shift = False
+
+        current_layer = 0
+
+        while not done:
+            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.update()
+
+                    if event.button == 5:
+                        v = -self.scale*np.cross(self.axis1,self.axis2)/10.0
+
+                        self.translate(v)
+
+                        self.update()
+
+                    if event.button == 1:
+                        clicked = True
+                        mouse_position = pygame.mouse.get_rel()
+
+                if event.type == MOUSEBUTTONUP:
+                    if event.button == 1:
+                        clicked = False
+
+                if event.type == MOUSEMOTION and clicked:
+                    movement = np.array(pygame.mouse.get_rel())
+
+                    if (movement == 0).all():
+                        continue
+
+                    length = np.linalg.norm(movement)
+
+                    mouse_direction = movement[0]*self.axis1 + movement[1]*self.axis2
+                    mouse_direction /= np.linalg.norm(mouse_direction)
+
+                    if shift:
+                        v = mouse_direction*self.scale*length/float(self.width)
+
+                        self.translate(v)
+
+                        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.update()
+
+                if event.type == KEYDOWN:
+                    if event.key == K_a:
+                        v = self.scale*self.axis1/10.0
+                        self.translate(v)
+                        self.update()
+
+                    if event.key == K_d:
+                        v = -self.scale*self.axis1/10.0
+                        self.translate(v)
+                        self.update()
+
+                    if event.key == K_w:
+                        v = self.scale*np.cross(self.axis1,self.axis2)/10.0
+                        self.translate(v)
+                        self.update()
+
+                    if event.key == K_s:
+                        v = -self.scale*np.cross(self.axis1,self.axis2)/10.0
+                        self.translate(v)
+                        self.update()
+
+                    if event.key == K_SPACE:
+                        v = self.scale*self.axis2/10.0
+                        self.translate(v)
+                        self.update()
+
+                    if event.key == K_LCTRL:
+                        v = -self.scale*self.axis2/10.0
+                        self.translate(v)
+                        self.update()
+
+                    if event.key == K_LSHIFT or event.key == K_RSHIFT:
+                        shift = True
+
+                    if event.key == K_ESCAPE:
+                        done = True
+                        break
+
+                    #if event.key == K_PAGEUP and load_bvh:
+                    #    try:
+                    #        if current_layer+1 >= len(bvhg):
+                    #            raise IndexError
+
+                    #        geometry = bvhg[current_layer+1]
+                    #        current_layer += 1
+
+                    #        geometry.load(module, color=True)
+                    #        update()
+                    #    except IndexError:
+                    #        print 'no further layers to view'
+
+                    #if event.key == K_PAGEDOWN and load_bvh:
+                    #    try:
+                    #        if current_layer-1 < 0:
+                    #            raise IndexError
+
+                    #        geometry = bvhg[current_layer-1]
+                    #        current_layer -= 1
+
+                    #        geometry.load(module, color=True)
+                    #        update()
+                    #    except IndexError:
+                    #        print 'no further layers to view'
+
+                    if event.key == K_F12:
+                        self.screenshot()
+
+                    if event.key == K_F5:
+                        self.render = not self.render
+                        self.update()
+
+                    if event.key == K_EQUALS:
+                        self.render_run_count += 1
+                        self.update()
+
+                    if event.key == K_MINUS:
+                        if self.render_run_count > 1:
+                            self.render_run_count -= 1
+                        self.update()
+
+                if event.type == KEYUP:
+                    if event.key == K_LSHIFT or event.key == K_RSHIFT:
+                        shift = False
+
+                if event.type == pygame.QUIT:
+                    done = True
+                    break
+
+        pygame.display.quit()
+
+if __name__ == '__main__':
+    import optparse
+    import inspect
+
+    import solids
+    import detectors
+    import scenes
+    from stl import mesh_from_stl
+    import src
+    from view import build
+
+    from pycuda.compiler import SourceModule
+
+    parser = optparse.OptionParser('%prog filename.stl')
+    parser.add_option('-b', '--bits', type='int', dest='bits',
+                      help='bits for z-ordering space axes', default=8)
+    #parser.add_option('-l', action='store_true', dest='load_bvh',
+    #                  help='load bounding volumes', default=False)
+    parser.add_option('-r', '--resolution', dest='resolution',
+                      help='specify resolution', default='800,600')
+    options, args = parser.parse_args()
+
+    if len(args) < 1:
+        sys.exit(parser.format_help())
+
+    size = [int(s) for s in options.resolution.split(',')]
+
+    if os.path.exists(args[0]):
+        root, ext = os.path.splitext(os.path.split(args[0])[1])
+
+        if ext.lower() == '.stl':
+            obj = mesh_from_stl(args[0])
+    else:
+        members = dict(inspect.getmembers(detectors) + inspect.getmembers(solids) + inspect.getmembers(scenes))
+
+        buildable_lookup = {}
+        for member in members.values():
+            if inspect.isfunction(member) and \
+                    hasattr(member, 'buildable') and member.buildable == True:
+                buildable_lookup[member.identifier] = member
+
+        if args[0] in buildable_lookup:
+            obj = buildable_lookup[args[0]]
+        else:
+            raise Exception("can't find object %s" % args[0])
+
+    from pycuda.tools import make_default_context
+
+    cuda.init()
+    context = make_default_context()
+
+    module = SourceModule(src.kernel, options=['-I' + src.dir], no_extern_c=True, cache_dir=False)
+    geometry = build(obj, options.bits)
+    geometry.load(module)
 
-        self.grid += (0,focal_length,0)
-        self.focal_point = np.zeros(3)
+    lock = Lock()
 
-    def position(self, position):
-        """Translate the camera to `position`."""
-        self.grid += position
-        self.focal_point += position
+    camera = Camera(geometry, module, context, lock, size)
 
-    def get_rays(self):
-        """Return the position and direction for each pixel ray."""
-        return self.grid, self.focal_point-self.grid
+    context.pop()
+    camera.start()
diff --git a/detectors/__init__.py b/detectors/__init__.py
index 2d7f250..18a1400 100644
--- a/detectors/__init__.py
+++ b/detectors/__init__.py
@@ -8,12 +8,12 @@ sys.path.append(dir + '/..')
 
 from view import buildable
 
-# from Water Cherenkov case study
-radius = 65.0/2
-height = 80.3
-nstrings = 229
-pmts_per_string = 90
-endcap_spacing = 0.8912
+# from LBNE document #3951
+radius = 63.30/2
+height = 76.60
+nstrings = 230
+pmts_per_string = 88
+endcap_spacing = 0.86
 
 @buildable('lbne')
 def build_lbne_200kton():
diff --git a/geometry.py b/geometry.py
index c76b089..1709a87 100644
--- a/geometry.py
+++ b/geometry.py
@@ -135,6 +135,9 @@ class Surface(object):
         else:
             value = np.tile(value, len(wavelengths))
 
+        if (np.asarray(value) < 0.0).any():
+            raise Exception('all probabilities must be >= 0.0')
+
         self.__dict__[name] = np.array(zip(wavelengths, value), dtype=np.float32)
 
 def interleave(arr, bits):
@@ -368,6 +371,9 @@ class Geometry(object):
         of the texture references.
         """
 
+        if not hasattr(self, 'mesh'):
+            raise RuntimeError('cannot load before a call to build().')
+
         set_wavelength_range = module.get_function('set_wavelength_range')
 
         set_wavelength_range(np.float32(standard_wavelengths[0]), np.float32(standard_wavelengths[-1]), np.float32(standard_wavelengths[1]-standard_wavelengths[0]), np.uint32(standard_wavelengths.size), block=(1,1,1), grid=(1,1))
diff --git a/optics.py b/optics.py
index 671d561..91f60cc 100644
--- a/optics.py
+++ b/optics.py
@@ -42,9 +42,10 @@ r7081hqe_photocathode.detect = \
               (710.0,  0.00)])
 # convert percent -> fraction
 r7081hqe_photocathode.detect[:,1] /= 100.0
-# photons not detected are absorbed
-#r7081hqe_photocathode.set('absorb', 1.0 - r7081hqe_photocathode.detect[:,1], wavelengths=r7081hqe_photocathode.detect[:,0])
-r7081hqe_photocathode.set('reflect_diffuse', 1.0 - r7081hqe_photocathode.detect[:,1], wavelengths=r7081hqe_photocathode.detect[:,0])
+# roughly the same amount of detected photons are absorbed without detection
+r7081hqe_photocathode.absorb = r7081hqe_photocathode.detect
+# remaining photons are diffusely reflected
+r7081hqe_photocathode.set('reflect_diffuse', 1.0 - r7081hqe_photocathode.detect[:,1] - r7081hqe_photocathode.absorb[:,1], wavelengths=r7081hqe_photocathode.detect[:,0])
 
 # water data comes from 'lightwater_sno' material in the SNO+ optics database
 water = Material('water')
diff --git a/render.py b/render.py
index 1890471..17a7289 100755
--- a/render.py
+++ b/render.py
@@ -21,6 +21,8 @@ def render(viewable, size=(800,600), name='', bits=8, make_movie=False):
 
     scale = np.linalg.norm(upper_bound-lower_bound)
 
+    from pycuda import autoinit
+
     print 'device %s' % autoinit.device.name()
 
     module = SourceModule(src.kernel, options=['-I' + src.dir], no_extern_c=True, cache_dir=False)
@@ -68,15 +70,15 @@ def render(viewable, size=(800,600), name='', bits=8, make_movie=False):
     rgb_lookup2[rgb_lookup2 > 1.0] = 1.0
     rgb_lookup2_gpu.set(rgb_lookup2.view(gpuarray.vec.float3))
 
-    camera = Camera(size)
+    #camera = Camera(size)
 
     point = np.array([0, diagonal*1.75, (lower_bound[2]+upper_bound[2])/2])
     axis1 = np.array([1,0,0], dtype=np.double)
     axis2 = np.array([0,0,1], dtype=np.double)
 
-    camera.position(point)
+    #camera.position(point)
 
-    origins, directions = camera.get_rays()
+    origins, directions = get_rays(point, size)
 
     origins_gpu = gpuarray.to_gpu(origins.astype(np.float32).view(gpuarray.vec.float3))
     directions_gpu = gpuarray.to_gpu(directions.astype(np.float32).view(gpuarray.vec.float3))
diff --git a/scenes/checkerboard.py b/scenes/checkerboard.py
index 1475a30..d2a9ec8 100644
--- a/scenes/checkerboard.py
+++ b/scenes/checkerboard.py
@@ -36,9 +36,9 @@ def build_checkerboard_scene(checkers_per_side=10, squares_per_checker=50):
 
     checkerboard = Solid(checkerboard_mesh, vacuum, vacuum, surface=checkerboard_surface, color=checkerboard_color)
 
-    sphere1 = Solid(sphere(theta=np.pi/128), water, vacuum)
-    sphere2 = Solid(sphere(theta=np.pi/128), vacuum, vacuum, surface=shiny_surface)
-    sphere3 = Solid(sphere(theta=np.pi/128), vacuum, vacuum, surface=lambertian_surface)
+    sphere1 = Solid(sphere(theta=np.pi/256), water, vacuum)
+    sphere2 = Solid(sphere(theta=np.pi/256), vacuum, vacuum, surface=shiny_surface)
+    sphere3 = Solid(sphere(theta=np.pi/256), vacuum, vacuum, surface=lambertian_surface)
 
     checkerboard_scene = Geometry()
     checkerboard_scene.add_solid(checkerboard, displacement=(0,0,-1.5))
diff --git a/solids/sno_cone.txt b/solids/sno_cone.txt
new file mode 100644
index 0000000..bf666cc
--- /dev/null
+++ b/solids/sno_cone.txt
@@ -0,0 +1,21 @@
+#DataThief /Users/stan/Downloads/./sno_pmt.png Tuesday 26-Jul-2011 3:33:52 PM
+-133.6402, 130.0879
+-133.673, 123.0691
+-133.7058, 116.0503
+-133.0342, 109.7357
+-132.3659, 102.7194
+-131.6944, 96.4048
+-130.3251, 89.3909
+-128.2481, 83.7831
+-126.8788, 76.7691
+-124.8018, 71.1613
+-122.7282, 64.8516
+-119.9503, 59.2462
+-117.1756, 52.9388
+-114.3976, 47.3335
+-111.6197, 41.7281
+-108.1407, 36.1251
+-104.6617, 30.5222
+-100.4816, 24.9216
+-97.5145, 22.2998
+
diff --git a/solids/sno_pmt.txt b/solids/sno_pmt.txt
new file mode 100644
index 0000000..b1fd494
--- /dev/null
+++ b/solids/sno_pmt.txt
@@ -0,0 +1,121 @@
+#DataThief /Users/stan/Downloads/./sno_pmt.png Tuesday 26-Jul-2011 3:31:47 PM
+-1.2088, -183.8939
+5.8016, -183.8698
+12.812, -183.8456
+19.8225, -183.8215
+26.8329, -183.7974
+33.8433, -183.7732
+38.0659, -180.2494
+38.0889, -175.3362
+36.7163, -169.0241
+36.749, -162.0053
+36.7818, -154.9865
+36.8145, -147.9677
+36.8473, -140.9489
+36.8801, -133.9301
+37.6139, -126.9089
+38.3477, -119.8877
+36.975, -113.5756
+35.6025, -107.2635
+35.6352, -100.2447
+34.9668, -93.2283
+34.2953, -86.9138
+35.7269, -80.592
+39.2616, -74.263
+43.4908, -69.3354
+48.4178, -65.1072
+54.0359, -62.9823
+60.3584, -60.1531
+62.491, -53.8289
+67.418, -49.6008
+72.3449, -45.3726
+77.2719, -41.1444
+82.9032, -36.212
+89.9136, -36.1878
+91.3418, -30.568
+94.8667, -26.3446
+97.7004, -20.0181
+99.8296, -14.3958
+101.2646, -7.3721
+101.9984, -0.3509
+101.3203, 4.5598
+100.652, 11.5762
+98.5783, 17.8859
+96.5013, 23.4937
+93.7267, 29.801
+90.2477, 35.404
+86.0644, 40.3026
+81.18, 45.1989
+76.2924, 49.3933
+71.4015, 52.8858
+65.8094, 56.376
+60.2176, 59.866
+54.6223, 62.6543
+48.326, 65.4401
+42.7275, 67.5264
+36.428, 69.6104
+30.1252, 70.9924
+23.8223, 72.3745
+16.8185, 73.7541
+10.5124, 74.4343
+3.5051, 75.1121
+-3.5051, 75.0879
+-10.5156, 75.0638
+-16.8316, 73.6383
+-23.8453, 72.9123
+-30.1612, 71.4868
+-36.4804, 69.3595
+-42.7995, 67.2322
+-48.4178, 65.1072
+-54.7403, 62.278
+-60.3617, 59.4512
+-65.983, 56.6244
+-71.6111, 52.3938
+-76.5348, 48.8675
+-82.1628, 44.6369
+-87.093, 39.7069
+-90.6212, 34.7816
+-94.1526, 29.1546
+-97.5145, 22.2998
+-99.1123, 17.9076
+-101.2449, 11.5834
+-101.9787, 4.5622
+-102.0115, -2.4566
+-101.3399, -8.7711
+-100.6716, -15.7875
+-98.5979, -22.0972
+-95.8167, -27.0007
+-93.0421, -33.308
+-88.8587, -38.2067
+-83.9744, -43.1029
+-79.09, -47.9992
+-74.1991, -51.4917
+-69.3082, -54.9843
+-63.7162, -58.4744
+-58.1243, -61.9645
+-52.5257, -64.0508
+-46.9338, -67.5409
+-42.0527, -73.1391
+-39.2715, -78.0426
+-36.4968, -84.3498
+-35.8252, -90.6644
+-37.2601, -97.688
+-37.2894, -104.0049
+-36.6213, -111.0213
+-40.1527, -116.6484
+-40.1822, -122.9654
+-40.2149, -129.9842
+-38.8423, -136.2963
+-38.8751, -143.3151
+-38.9078, -150.3338
+-38.9405, -157.3527
+-38.9734, -164.3715
+-39.0061, -171.3902
+-41.8365, -177.015
+-41.1649, -183.3295
+-34.8588, -184.0097
+-27.8484, -183.9855
+-20.838, -183.9614
+-13.8276, -183.9373
+-6.8171, -183.9132
+-3.3119, -183.9011
diff --git a/src/kernel.cu b/src/kernel.cu
index 0dc5639..ffd2d88 100644
--- a/src/kernel.cu
+++ b/src/kernel.cu
@@ -10,7 +10,7 @@
 #include "alpha.h"
 
 #define RED_WAVELENGTH 685
-#define BLUE_WAVELENGTH 465
+#define BLUE_WAVELENGTH 445
 #define GREEN_WAVELENGTH 545
 
 __device__ void fAtomicAdd(float *addr, float data)
@@ -147,7 +147,7 @@ __global__ void build_rgb_lookup(int nthreads, curandState *rng_states, float3 *
 		}
 
 		p = seed;
-		p.wavelength = BLUE_WAVELENGTH;
+		p.wavelength = GREEN_WAVELENGTH;
 
 		to_diffuse(p, s, rng, max_steps);
 
@@ -164,7 +164,7 @@ __global__ void build_rgb_lookup(int nthreads, curandState *rng_states, float3 *
 		}
 
 		p = seed;
-		p.wavelength = GREEN_WAVELENGTH;
+		p.wavelength = BLUE_WAVELENGTH;
 
 		to_diffuse(p, s, rng, max_steps);
 
@@ -227,7 +227,7 @@ __global__ void render(int nthreads, curandState *rng_states, float3 *positions,
 		}
 
 		p = seed;
-		p.wavelength = BLUE_WAVELENGTH;
+		p.wavelength = GREEN_WAVELENGTH;
 
 		to_diffuse(p, s, rng, max_steps);
 
@@ -244,7 +244,7 @@ __global__ void render(int nthreads, curandState *rng_states, float3 *positions,
 		}
 
 		p = seed;
-		p.wavelength = GREEN_WAVELENGTH;
+		p.wavelength = BLUE_WAVELENGTH;
 
 		to_diffuse(p, s, rng, max_steps);
 
diff --git a/view.py b/view.py
index 11c4b57..96f7530 100755
--- a/view.py
+++ b/view.py
@@ -9,12 +9,12 @@ import pygame
 from pygame.locals import *
 
 import src
-from camera import Camera
+from camera import get_rays
 from geometry import Mesh, Solid, Geometry
 from transform import rotate
 from optics import *
 
-from pycuda import autoinit
+#from pycuda import autoinit
 from pycuda.compiler import SourceModule
 from pycuda import gpuarray
 import pycuda.driver as cuda
@@ -158,6 +158,8 @@ def view(viewable, size=(800,600), name='', bits=8, load_bvh=False):
 
     scale = np.linalg.norm(upper_bound-lower_bound)
 
+    from pycuda import autoinit
+
     print 'device %s' % autoinit.device.name()
 
     module = SourceModule(src.kernel, options=['-I' + src.dir], no_extern_c=True, cache_dir=False)
@@ -171,7 +173,7 @@ def view(viewable, size=(800,600), name='', bits=8, load_bvh=False):
     screen = pygame.display.set_mode(size)
     pygame.display.set_caption(name)
 
-    camera = Camera(size)
+    #camera = Camera(size)
 
     diagonal = np.linalg.norm(upper_bound-lower_bound)
 
@@ -179,9 +181,9 @@ def view(viewable, size=(800,600), name='', bits=8, load_bvh=False):
     axis1 = np.array([1,0,0], dtype=np.double)
     axis2 = np.array([0,0,1], dtype=np.double)
 
-    camera.position(point)
+    #camera.position(point)
 
-    origins, directions = camera.get_rays()
+    origins, directions = get_rays(point, size)
 
     origins_gpu = gpuarray.to_gpu(origins.astype(np.float32).view(gpuarray.vec.float3))