diff options
| -rwxr-xr-x | camera.py | 476 | ||||
| -rw-r--r-- | detectors/__init__.py | 4 | ||||
| -rw-r--r-- | detectors/lbne.py | 11 | ||||
| -rw-r--r-- | geometry.py | 2 | ||||
| -rw-r--r-- | gpu.py | 32 | ||||
| -rwxr-xr-x | sim.py | 2 | ||||
| -rw-r--r-- | solids/__init__.py | 6 | ||||
| -rw-r--r-- | solids/pmts.py | 22 | ||||
| -rw-r--r-- | src/mesh.h | 17 | ||||
| -rw-r--r-- | tools.py | 10 | ||||
| -rwxr-xr-x | view.py | 20 |
11 files changed, 311 insertions, 291 deletions
@@ -2,12 +2,12 @@ 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 color import map_wavelength, map_to_color +from gpu import GPU, CUDAFuncs +from tools import timeit from transform import rotate import pygame @@ -21,13 +21,9 @@ 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 +from histogram import * + +import matplotlib.pyplot as plt def encode_movie(dir): root, ext = 'movie', 'avi' @@ -72,24 +68,24 @@ def get_rays(position, size = (800, 600), film_size = (0.035, 0.024), focal_leng return grid, focal_point-grid class Camera(Thread): - def __init__(self, geometry, module, context, lock=None, size=(800,600), - spnav=False): + def __init__(self, geometry, size=(800,600), device_id=None, spnav=False): Thread.__init__(self) self.geometry = geometry - self.module = module - self.context = context + self.device_id = device_id + self.size = size + self.spnav = spnav if spnav: import spnav as spnav_module self.spnav_module = spnav_module - if lock is None: - self.lock = Lock() - else: - self.lock = lock + def init_gpu(self): + self.gpu = GPU(self.device_id) + self.gpu.load_geometry(geometry) + + self.kernels = CUDAFuncs(self.gpu.module, ['ray_trace', 'rotate', 'rotate_around_point', 'translate', 'update_xyz_lookup', 'update_xyz_image', 'process_image', 'init_rng']) - self.size = size self.width, self.height = size pygame.init() @@ -97,18 +93,6 @@ class Camera(Thread): pygame.display.set_caption('') self.clock = pygame.time.Clock() - 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.rotate_around_point_kernel = self.module.get_function('rotate_around_point') - self.translate_kernel = self.module.get_function('translate') - 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() - lower_bound, upper_bound = self.geometry.mesh.get_bounds() self.scale = np.linalg.norm(upper_bound-lower_bound) @@ -120,12 +104,9 @@ class Camera(Thread): origins, directions = get_rays(self.point, self.size) - 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.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.movie = False self.movie_index = 0 @@ -134,16 +115,12 @@ class Camera(Thread): @timeit def initialize_render(self): - with self.lock: - self.context.push() - print 'beginning initialize_render()' - 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() - self.context.pop() + self.rng_states_gpu = cuda.mem_alloc(self.width*self.height*sizeof('curandStateXORWOW', '#include <curand_kernel.h>')) + self.kernels.init_rng(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() self.source_position = self.point @@ -152,54 +129,39 @@ class Camera(Thread): self.max_steps = 10 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() + 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.nlookup_calls = 0 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)) + for i in range(self.xyz_lookup1_gpu.size//(self.width*self.height)+1): + self.kernels.update_xyz_lookup(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)) - 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.kernels.update_xyz_lookup(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() + for i in range(self.xyz_lookup1_gpu.size//(self.width*self.height)+1): + self.kernels.update_xyz_lookup(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.nlookup_calls += 1 def clear_image(self): - with self.lock: - self.context.push() - self.image_gpu.fill(gpuarray.vec.make_float3(0.0,0.0,0.0)) - self.context.pop() + self.image_gpu.fill(gpuarray.vec.make_float3(0.0,0.0,0.0)) self.nimages = 0 def update_image(self): - with self.lock: - self.context.push() - 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.kernels.update_xyz_image(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)) + self.kernels.update_xyz_image(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)) - 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() + self.kernels.update_xyz_image(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.nimages += 1 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() + self.kernels.process_image(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)) def screenshot(self, dir='', start=0): root, ext = 'screenshot', 'png' @@ -214,15 +176,12 @@ class Camera(Thread): 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, 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.kernels.rotate(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.kernels.rotate(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() + self.point = rotate(self.point, phi, n) + self.axis1 = rotate(self.axis1, phi, n) + self.axis2 = rotate(self.axis2, phi, n) if self.render: self.clear_image() @@ -230,11 +189,8 @@ class Camera(Thread): self.update() def rotate_around_point(self, phi, n, point, redraw=True): - 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.kernels.rotate_around_point(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.kernels.rotate(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.axis1 = rotate(self.axis1, phi, n) self.axis2 = rotate(self.axis2, phi, n) @@ -246,12 +202,9 @@ class Camera(Thread): self.update() def translate(self, v, redraw=True): - 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.kernels.translate(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() + self.point += v if redraw: if self.render: @@ -266,38 +219,135 @@ class Camera(Thread): 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() + self.kernels.ray_trace(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)) - with self.lock: - self.context.push() - pygame.surfarray.blit_array(self.screen, self.pixels_gpu.get().reshape(self.size)) - pygame.display.flip() - self.context.pop() + pygame.surfarray.blit_array(self.screen, self.pixels_gpu.get().reshape(self.size)) + pygame.display.flip() if self.movie: self.screenshot(self.movie_dir, self.movie_index) self.movie_index += 1 + def process_event(self, event): + if event.type == MOUSEBUTTONDOWN: + if event.button == 4: + v = self.scale*np.cross(self.axis1,self.axis2)/10.0 + self.translate(v) + + elif event.button == 5: + v = -self.scale*np.cross(self.axis1,self.axis2)/10.0 + self.translate(v) + + elif event.button == 1: + mouse_position = pygame.mouse.get_rel() + self.clicked = True + + elif event.type == MOUSEBUTTONUP: + if event.button == 1: + self.clicked = False + + elif event.type == MOUSEMOTION and self.clicked: + movement = np.array(pygame.mouse.get_rel()) + + if (movement == 0).all(): + return + + length = np.linalg.norm(movement) + + mouse_direction = movement[0]*self.axis1 + movement[1]*self.axis2 + mouse_direction /= np.linalg.norm(mouse_direction) + + if pygame.key.get_mods() & (KMOD_LSHIFT | KMOD_RSHIFT): + v = mouse_direction*self.scale*length/float(self.width) + self.translate(v) + else: + phi = np.float32(2*np.pi*length/float(self.width)) + n = rotate(mouse_direction, np.pi/2, -np.cross(self.axis1,self.axis2)) + + if pygame.key.get_mods() & KMOD_LCTRL: + self.rotate_around_point(phi, n, self.point) + else: + self.rotate(phi, n) + + elif event.type == KEYDOWN: + if event.key == K_a: + v = self.scale*self.axis1/10.0 + self.translate(v) + + elif event.key == K_d: + v = -self.scale*self.axis1/10.0 + self.translate(v) + + elif event.key == K_w: + v = self.scale*np.cross(self.axis1,self.axis2)/10.0 + self.translate(v) + + elif event.key == K_s: + v = -self.scale*np.cross(self.axis1,self.axis2)/10.0 + self.translate(v) + + elif event.key == K_SPACE: + v = self.scale*self.axis2/10.0 + self.translate(v) + + elif event.key == K_F6: + self.clear_xyz_lookup() + self.clear_image() + self.source_position = self.point + + elif event.key == K_p: + for i in range(100): + self.update_xyz_lookup(self.point) + self.source_position = self.point + + elif event.key == K_ESCAPE: + self.done = True + return + + elif event.key == K_F12: + self.screenshot() + + elif event.key == K_F5: + if not hasattr(self, 'rng_states_gpu'): + self.initialize_render() + + self.render = not self.render + self.clear_image() + self.update() + + elif 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 + + elif event.type == pygame.QUIT: + self.done = True + return + def run(self): + self.init_gpu() + if self.spnav: self.spnav_module.spnav_open() self.update() - done = False - clicked = False - shift = False - ctrl = False + self.done = False + self.clicked = False + #shift = False + #ctrl = False #current_layer = 0 - while not done: + while not self.done: self.clock.tick(20) - if self.render and not clicked and not pygame.event.peek(KEYDOWN): + if self.render and not self.clicked and not pygame.event.peek(KEYDOWN): self.update() # Space Navigator controls @@ -353,156 +403,70 @@ class Camera(Thread): 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) - - if event.button == 5: - v = -self.scale*np.cross(self.axis1,self.axis2)/10.0 - self.translate(v) - - 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 + self.process_event(event) - 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) - else: - phi = np.float32(2*np.pi*length/float(self.width)) - n = rotate(mouse_direction, np.pi/2, -np.cross(self.axis1,self.axis2)) - - 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) - - if event.key == K_d: - v = -self.scale*self.axis1/10.0 - self.translate(v) - - if event.key == K_w: - v = self.scale*np.cross(self.axis1,self.axis2)/10.0 - self.translate(v) - - if event.key == K_s: - v = -self.scale*np.cross(self.axis1,self.axis2)/10.0 - self.translate(v) - - if event.key == K_SPACE: - 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.source_position = self.point - - if event.key == K_p: - for i in range(100): - 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 - - #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 + if self.movie: + encode_movie(self.movie_dir) - # geometry.load(module, color=True) - # update() - # except IndexError: - # print 'no further layers to view' + pygame.display.quit() + if self.spnav: + self.spnav_module.spnav_close() - if event.key == K_F12: - self.screenshot() + del self.gpu - if event.key == K_F5: - if not hasattr(self, 'rng_states_gpu'): - self.initialize_render() +class EventViewer(Camera): + def __init__(self, geometry, filename, **kwargs): + Camera.__init__(self, geometry, **kwargs) - self.render = not self.render - self.clear_image() - self.update() + import ROOT - 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 + self.f = ROOT.TFile(filename) + self.T = self.f.Get('T') + self.T.GetEntry(0) - if event.type == KEYUP: - if event.key == K_LSHIFT or event.key == K_RSHIFT: - shift = False + @timeit + def color_hit_pmts(self): + self.gpu.reset_colors() - if event.key == K_LCTRL or event.key == K_RCTRL: - ctrl = False + solid_ids = np.empty(len(self.T.ev.channel), np.uint32) + t = np.empty(len(self.T.ev.channel), np.float32) + q = np.empty(len(self.T.ev.channel), np.float32) - if event.type == pygame.QUIT: - done = True - break + for i, channel in enumerate(self.T.ev.channel): + solid_ids[i] = channel.channel_id + t[i] = channel.t + q[i] = channel.q - if self.movie: - encode_movie(self.movie_dir) + self.gpu.color_solids(solid_ids, map_to_color(t)) + self.update() - pygame.display.quit() - if self.spnav: - self.spnav_module.spnav_close() + plt.clf() + plt.hist(t*1e9, 100) + plt.xlabel('Time (ns)') + plt.ylabel('Entries') + plt.draw() + + def process_event(self, event): + if event.type == KEYDOWN: + if event.key == K_PAGEUP: + entry = self.T.GetReadEntry() + if entry < self.T.GetEntries() - 1: + self.T.GetEntry(entry+1) + self.color_hit_pmts() + return + + if event.key == K_PAGEDOWN: + entry = self.T.GetReadEntry() + if entry > 0: + self.T.GetEntry(entry-1) + self.color_hit_pmts() + return + + Camera.process_event(self, event) + def run(self): + plt.ion() + Camera.run(self) if __name__ == '__main__': import optparse @@ -512,11 +476,8 @@ if __name__ == '__main__': 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=10) @@ -524,6 +485,7 @@ if __name__ == '__main__': # help='load bounding volumes', default=False) parser.add_option('-r', '--resolution', dest='resolution', help='specify resolution', default='800,600') + parser.add_option('-i', dest='io_file', default=None) parser.add_option('--spnav', action='store_true', dest='spnav', help='activate Space Navigator support', default=False) options, args = parser.parse_args() @@ -538,6 +500,7 @@ if __name__ == '__main__': if ext.lower() in ('.stl', '.bz2'): obj = mesh_from_stl(args[0]) + else: members = dict(inspect.getmembers(detectors) + inspect.getmembers(solids) + inspect.getmembers(scenes)) @@ -552,16 +515,9 @@ if __name__ == '__main__': else: raise Exception("can't find object %s" % args[0]) - from pycuda.tools import make_default_context - from gpu import * - geometry = build(obj, options.bits) - - lock = Lock() - - gpu = GPU() - gpu.load_geometry(geometry) - gpu.context.pop() - camera = Camera(geometry, gpu.module, gpu.context, lock, size=size, - spnav=options.spnav) + if options.io_file is not None: + camera = EventViewer(geometry, options.io_file, size=size, spnav=options.spnav) + else: + camera = Camera(geometry, size, spnav=options.spnav) camera.start() diff --git a/detectors/__init__.py b/detectors/__init__.py index 043bdbd..8550751 100644 --- a/detectors/__init__.py +++ b/detectors/__init__.py @@ -17,6 +17,10 @@ nstrings = 230 pmts_per_string = 88 endcap_spacing = 0.86 +@buildable('lbne_event_view') +def build_lbne_200kton_event_view(): + return build_lbne(radius, height, nstrings, pmts_per_string, endcap_spacing, physical_model=False) + @buildable('lbne') def build_lbne_200kton(): return build_lbne(radius, height, nstrings, pmts_per_string, endcap_spacing) diff --git a/detectors/lbne.py b/detectors/lbne.py index 8fe2c18..3e8ae7b 100644 --- a/detectors/lbne.py +++ b/detectors/lbne.py @@ -6,21 +6,24 @@ dir = os.path.split(os.path.realpath(__file__))[0] sys.path.append(dir + '/..') from geometry import * -from solids import build_12inch_pmt, build_12inch_pmt_with_lc +from solids import build_12inch_pmt, build_12inch_pmt_with_lc, build_12inch_pmt_shell from optics import * from transform import rotate, make_rotation_matrix from itertools import product from make import cylinder -def build_lbne(radius, height, nstrings, pmts_per_string, endcap_spacing): - pmt = build_12inch_pmt_with_lc() +def build_lbne(radius, height, nstrings, pmts_per_string, endcap_spacing, physical_model=True): + if physical_model: + pmt = build_12inch_pmt_with_lc() + else: + pmt = build_12inch_pmt_shell() lbne = Geometry() # outer cylinder cylinder_mesh = cylinder(radius, radius, height+height/(pmts_per_string-1), theta=(2*np.pi/nstrings)/4) cylinder_mesh.vertices = rotate(cylinder_mesh.vertices, np.pi/2, (-1,0,0)) - lbne.add_solid(Solid(cylinder_mesh, water, vacuum, black_surface, 0x990000ff)) + lbne.add_solid(Solid(cylinder_mesh, water, vacuum, black_surface, 0xff0000ff)) lbne.pmtids = [] diff --git a/geometry.py b/geometry.py index 1b85c51..37367e5 100644 --- a/geometry.py +++ b/geometry.py @@ -4,7 +4,7 @@ import pycuda.driver as cuda from pycuda import gpuarray from itertools import * from itertoolset import * -from camera import timeit +from tools import timeit # all material/surface properties are interpolated at these # wavelengths when they are sent to the gpu @@ -7,6 +7,7 @@ import pycuda.driver as cuda from pycuda import gpuarray from copy import copy from itertools import izip +from tools import timeit import src from geometry import standard_wavelengths @@ -42,18 +43,16 @@ class GPU(object): '''Initialize a GPU context on the specified device. If device_id==None, the default device is used.''' - if device_id == None: + if device_id is None: self.context = make_default_context() else: - device = pycuda.driver.Device(device_id) + print 'device_id = %s' % device_id + device = cuda.Device(device_id) self.context = device.make_context() print 'device %s' % self.context.get_device().name() self.module = SourceModule(src.kernel, options=['-I' + src.dir], no_extern_c=True, cache_dir=False) - self.geo_funcs = CUDAFuncs(self.module, ['set_wavelength_range', 'set_material', - 'set_surface', - 'set_global_mesh_variables']) - - + self.geo_funcs = CUDAFuncs(self.module, ['set_wavelength_range', 'set_material', 'set_surface', 'set_global_mesh_variables', 'color_solids']) + self.prop_funcs = CUDAFuncs(self.module, ['init_rng', 'propagate']) self.nblocks = 64 self.max_nthreads = 200000 @@ -154,16 +153,21 @@ class GPU(object): self.node_map_tex.set_format(cuda.array_format.UNSIGNED_INT32, 1) self.node_length_tex.set_format(cuda.array_format.UNSIGNED_INT32, 1) + self.geometry = geometry + self.print_device_usage() - def color_solids(self, solid_ids, colors): - solid_id_map = self.solid_id_map_gpu.get() - triangle_colors = self.colors_gpu.get() + @timeit + def reset_colors(self): + self.colors_gpu.set(self.geometry.colors.astype(np.uint32)) - for i, color in izip(solid_ids, colors): - triangle_colors[solid_id_map == i] = color + @timeit + def color_solids(self, solid_ids, colors): + solid_ids_gpu = gpuarray.to_gpu(np.array(solid_ids, dtype=np.int32)) + solid_colors_gpu = gpuarray.to_gpu(np.array(colors, dtype=np.uint32)) - self.colors_gpu.set(triangle_colors) + self.geo_funcs.color_solids(np.int32(solid_ids_gpu.size), np.uint32(self.triangles_gpu.size), self.solid_id_map_gpu, solid_ids_gpu, solid_colors_gpu, block=(self.nblocks,1,1), grid=(solid_ids_gpu.size//self.nblocks+1,1)) + self.context.synchronize() def setup_propagate(self, seed=1): self.rng_states_gpu = cuda.mem_alloc(self.max_nthreads*sizeof('curandStateXORWOW', '#include <curand_kernel.h>')) @@ -277,5 +281,5 @@ class GPU(object): def get_hits(self): return self.earliest_time_gpu.get() - def shutdown(self): + def __del__(self): self.context.pop() @@ -136,5 +136,3 @@ if __name__ == '__main__': f.Close() print >>sys.stderr, 'Done. %1.1f events/sec, %1.0f photons/sec.' % (options.nevents/(end_sim - start_sim), nphotons/(end_sim - start_sim)) - - gpu_worker.shutdown() diff --git a/solids/__init__.py b/solids/__init__.py index 479bacd..46f3119 100644 --- a/solids/__init__.py +++ b/solids/__init__.py @@ -1,6 +1,6 @@ import numpy as np -from pmts import build_pmt, build_light_collector, build_light_collector_from_file +from pmts import build_pmt, build_light_collector, build_light_collector_from_file, build_pmt_shell import os import sys @@ -15,6 +15,10 @@ from view import buildable def build_12inch_pmt(outer_material=water, theta=np.pi/8): return build_pmt(dir + '/hamamatsu_12inch.txt', 0.003, outer_material, theta) +@buildable('12inch_pmt_shell') +def build_12inch_pmt_shell(outer_material=water, theta=np.pi/8): + return build_pmt_shell(dir + '/hamamatsu_12inch.txt') + # from Jelena Maricic lc_12inch_a = 0.16597 lc_12inch_b = 0.584525 diff --git a/solids/pmts.py b/solids/pmts.py index 22d9d67..a30d007 100644 --- a/solids/pmts.py +++ b/solids/pmts.py @@ -27,7 +27,24 @@ def build_light_collector(pmt, a, b, d, rmin, rmax, npoints=10): lc_mesh = rotate_extrude(lc_radii, lc_profile + lc_offset, pmt.theta) - return Solid(lc_mesh, pmt.outer_material, pmt.outer_material, color=0xff0000, surface=shiny_surface) + return Solid(lc_mesh, pmt.outer_material, pmt.outer_material, surface=shiny_surface) + +def build_pmt_shell(filename, outer_material=water, theta=np.pi/8): + profile = read_csv(filename) + + # slice profile in half + profile = profile[profile[:,0] < 0] + profile[:,0] = -profile[:,0] + # order profile from base to face + profile = profile[np.argsort(profile[:,1])] + # set x coordinate to 0.0 for first and last profile along the profile + # so that the mesh is closed + profile[0,0] = 0.0 + profile[-1,0] = 0.0 + # convert mm -> m + profile /= 1000.0 + + return Solid(rotate_extrude(profile[:,0], profile[:,1], theta), glass, outer_material, color=0xeeffffff) def build_pmt(filename, glass_thickness, outer_material=water, theta=np.pi/8): profile = read_csv(filename) @@ -72,6 +89,5 @@ def build_light_collector_from_file(filename, outer_material, theta=np.pi/24): profile /= 1000.0 mesh = rotate_extrude(profile[:,0], profile[:,1], theta) - solid = Solid(mesh, outer_material, outer_material, color=0xFF0000, - surface=shiny_surface) + solid = Solid(mesh, outer_material, outer_material, surface=shiny_surface) return solid @@ -148,6 +148,23 @@ __global__ void set_colors(unsigned int *colors) g_colors = colors; } +__global__ void color_solids(int nthreads, unsigned int ntriangles, int *solid_id_map, int *solid_ids, unsigned int *solid_colors) +{ + int id = blockIdx.x*blockDim.x + threadIdx.x; + + if (id >= nthreads) + return; + + int solid_id = solid_ids[id]; + unsigned int color = solid_colors[id]; + + for (int i=0; i < ntriangles; i++) + { + if (solid_id_map[i] == solid_id) + g_colors[i] = color; + } +} + } // extern "c" #endif @@ -1,4 +1,14 @@ import numpy as np +import time +import datetime + +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 read_csv(filename): """Return an array of comma-separated values from `filename`.""" @@ -13,6 +13,7 @@ from camera import get_rays from geometry import Mesh, Solid, Geometry from transform import rotate from optics import * +from gpu import * #from pycuda import autoinit from pycuda.compiler import SourceModule @@ -142,6 +143,8 @@ def view(viewable, size=(800,600), name='', bits=8, load_bvh=False): - move: shift+click and drag the mouse """ + gpu = GPU() + geometry = build(viewable, bits) if load_bvh: @@ -158,12 +161,15 @@ def view(viewable, size=(800,600), name='', bits=8, load_bvh=False): scale = np.linalg.norm(upper_bound-lower_bound) - from pycuda import autoinit + #from pycuda import autoinit + + #print 'device %s' % autoinit.device.name() - print 'device %s' % autoinit.device.name() + module = gpu.module - module = SourceModule(src.kernel, options=['-I' + src.dir], no_extern_c=True, cache_dir=False) - geometry.load(module) + #module = SourceModule(src.kernel, options=['-I' + src.dir], no_extern_c=True)#, cache_dir=False) + #geometry.load(module) + gpu.load_geometry(geometry) cuda_raytrace = module.get_function('ray_trace') cuda_rotate = module.get_function('rotate') cuda_translate = module.get_function('translate') @@ -292,7 +298,8 @@ def view(viewable, size=(800,600), name='', bits=8, load_bvh=False): geometry = bvhg[current_layer+1] current_layer += 1 - geometry.load(module, color=True) + gpu.load_geometry(geometry) + #geometry.load(module, color=True) update() except IndexError: print 'no further layers to view' @@ -305,7 +312,8 @@ def view(viewable, size=(800,600), name='', bits=8, load_bvh=False): geometry = bvhg[current_layer-1] current_layer -= 1 - geometry.load(module, color=True) + gpu.load_geometry(geometry) + #geometry.load(module, color=True) update() except IndexError: print 'no further layers to view' |