summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--geometry.py49
-rw-r--r--gputhread.py27
-rw-r--r--itertoolset.py8
-rw-r--r--make.py6
-rw-r--r--materials/__init__.py29
-rw-r--r--models/portal_gun/Portal_Gun_Barrel_Insert.STLbin0 -> 122284 bytes
-rw-r--r--models/portal_gun/Portal_Gun_Bottom_Cover.STLbin0 -> 3691384 bytes
-rw-r--r--models/portal_gun/Portal_Gun_Center_Tube.STLbin0 -> 45884 bytes
-rw-r--r--models/portal_gun/Portal_Gun_Claw.STLbin0 -> 248084 bytes
-rw-r--r--models/portal_gun/Portal_Gun_Left_Half.STLbin0 -> 2139984 bytes
-rw-r--r--models/portal_gun/Portal_Gun_Right_Half.STLbin0 -> 2387784 bytes
-rw-r--r--models/portal_gun/Portal_Gun_Top_Claw_Mount.STLbin0 -> 202884 bytes
-rw-r--r--optics.py148
-rw-r--r--profiles/hamamatsu_12inch.py112
-rw-r--r--profiles/hamamatsu_12inch.txt55
-rw-r--r--profiles/test.py3
-rwxr-xr-xrender.py468
-rw-r--r--scenes/checkerboard.py58
-rw-r--r--solids/r11708.py2
-rw-r--r--solids/r11708_cut.py2
-rw-r--r--src/kernel.cu450
-rw-r--r--src/materials.h16
-rw-r--r--threadtest.py2
-rwxr-xr-xview.py194
24 files changed, 1528 insertions, 101 deletions
diff --git a/geometry.py b/geometry.py
index f0f5c20..08ce291 100644
--- a/geometry.py
+++ b/geometry.py
@@ -1,3 +1,4 @@
+import sys
import numpy as np
import pycuda.driver as cuda
from pycuda import gpuarray
@@ -119,9 +120,15 @@ class Surface(object):
def __init__(self, name='none'):
self.name = name
- self.absorption = None
- self.reflection_diffuse = None
- self.reflection_specular = None
+ self.set('detect', 0)
+ self.set('absorb', 0)
+ self.set('reflect_diffuse', 0)
+ self.set('reflect_specular', 0)
+
+ #self.detect = None
+ #self.absorb = None
+ #self.reflect_diffuse = None
+ #self.reflect_specular = None
def set(self, name, value, wavelengths=standard_wavelengths):
if np.iterable(value):
@@ -202,7 +209,7 @@ class Geometry(object):
return len(self.solids)-1
- def build(self, bits=8):
+ def build(self, bits=8, shift=3):
offsets = [ (0,0) ]
for solid in self.solids:
offsets.append( (offsets[-1][0] + len(solid.mesh.vertices),
@@ -246,7 +253,10 @@ class Geometry(object):
self.surface_index = \
np.fromiter(imap(surface_lookup.get, chain(*[solid.surface for solid in self.solids])), dtype=np.int32)[reorder]
- self.surface_index[self.surface_index == self.unique_surfaces.index(None)] = -1
+ try:
+ self.surface_index[self.surface_index == self.unique_surfaces.index(None)] = -1
+ except ValueError:
+ pass
self.colors = \
np.concatenate([solid.color for solid in self.solids])[reorder]
@@ -262,6 +272,9 @@ class Geometry(object):
self.upper_bounds = np.empty((unique_zvalues.size,3), dtype=np.float32)
self.node_map = np.empty(unique_zvalues.size, dtype=np.uint32)
self.node_length = np.empty(unique_zvalues.size, dtype=np.uint32)
+ self.layers = np.empty(unique_zvalues.size, dtype=np.uint32)
+
+ layer = 0
for i, z in enumerate(unique_zvalues):
i1 = np.searchsorted(zvalues_mesh, z)
@@ -277,6 +290,7 @@ class Geometry(object):
self.node_map[i] = i1
self.node_length[i] = i2-i1
+ self.layers[i] = layer
zvalues[i] = z
@@ -284,8 +298,10 @@ class Geometry(object):
begin_last_layer = 0
+ layer += 1
+
while True:
- bit_shifted_zvalues = zvalues >> 1
+ bit_shifted_zvalues = zvalues >> shift
unique_bit_shifted_zvalues = np.unique(bit_shifted_zvalues)
zvalues = np.empty(unique_bit_shifted_zvalues.size, dtype=np.uint64)
@@ -299,6 +315,8 @@ class Geometry(object):
self.node_length.resize(\
self.node_length.size+unique_bit_shifted_zvalues.size)
+ self.layers.resize(self.layers.size+unique_bit_shifted_zvalues.size)
+
for i, z in enumerate(unique_bit_shifted_zvalues):
i1 = np.searchsorted(bit_shifted_zvalues, z) + \
begin_last_layer
@@ -321,9 +339,12 @@ class Geometry(object):
self.node_map[i] = i1
self.node_length[i] = i2-i1
+ self.layers[i] = layer
begin_last_layer += bit_shifted_zvalues.size
+ layer += 1
+
if unique_bit_shifted_zvalues.size == 1:
break
@@ -360,15 +381,17 @@ class Geometry(object):
if surface is None:
continue
- absorption = np.interp(standard_wavelengths, surface.absorption[:,0], surface.absorption[:,1]).astype(np.float32)
- reflection_diffuse = np.interp(standard_wavelengths, surface.reflection_diffuse[:,0], surface.reflection_diffuse[:,1]).astype(np.float32)
- reflection_specular = np.interp(standard_wavelengths, surface.reflection_specular[:,0], surface.reflection_specular[:,1]).astype(np.float32)
+ detect = np.interp(standard_wavelengths, surface.detect[:,0], surface.detect[:,1]).astype(np.float32)
+ absorb = np.interp(standard_wavelengths, surface.absorb[:,0], surface.absorb[:,1]).astype(np.float32)
+ reflect_diffuse = np.interp(standard_wavelengths, surface.reflect_diffuse[:,0], surface.reflect_diffuse[:,1]).astype(np.float32)
+ reflect_specular = np.interp(standard_wavelengths, surface.reflect_specular[:,0], surface.reflect_specular[:,1]).astype(np.float32)
- surface.absorption_gpu = cuda.to_device(absorption)
- surface.reflection_diffuse_gpu = cuda.to_device(reflection_diffuse)
- surface.reflection_specular_gpu = cuda.to_device(reflection_specular)
+ surface.detect_gpu = cuda.to_device(detect)
+ surface.absorb_gpu = cuda.to_device(absorb)
+ surface.reflect_diffuse_gpu = cuda.to_device(reflect_diffuse)
+ surface.reflect_specular_gpu = cuda.to_device(reflect_specular)
- set_surface(np.int32(i), surface.absorption_gpu, surface.reflection_diffuse_gpu, surface.reflection_specular_gpu, block=(1,1,1), grid=(1,1))
+ set_surface(np.int32(i), surface.detect_gpu, surface.absorb_gpu, surface.reflect_diffuse_gpu, surface.reflect_specular_gpu, block=(1,1,1), grid=(1,1))
vertices = np.empty(len(self.mesh.vertices), dtype=gpuarray.vec.float3)
vertices['x'] = self.mesh.vertices[:,0]
diff --git a/gputhread.py b/gputhread.py
index 06921d5..f9e6dae 100644
--- a/gputhread.py
+++ b/gputhread.py
@@ -33,8 +33,6 @@ class GPUThread(threading.Thread):
module = SourceModule(src.kernel, options=['-I' + src.dir], no_extern_c=True, cache_dir=False)
propagate = module.get_function('propagate')
- fill_float = module.get_function('fill_float')
- fill_float3 = module.get_function('fill_float3')
fill_uniform = module.get_function('fill_uniform')
fill_uniform_sphere = module.get_function('fill_uniform_sphere')
@@ -42,7 +40,7 @@ class GPUThread(threading.Thread):
rng_states_gpu = cuda.mem_alloc(self.max_nthreads*sizeof('curandStateXORWOW', '#include <curand_kernel.h>'))
init_rng(np.int32(self.max_nthreads), rng_states_gpu, np.int32(self.device_id), np.int32(0), block=(self.nblocks,1,1), grid=(self.max_nthreads//self.nblocks+1,1))
- texrefs = self.geometry.load(module)
+ self.geometry.load(module)
daq_module = SourceModule(src.daq, options=['-I' + src.dir], no_extern_c=True, cache_dir=False)
reset_earliest_time_int = daq_module.get_function('reset_earliest_time_int')
@@ -64,24 +62,23 @@ class GPUThread(threading.Thread):
gpu_kwargs = {'block': (self.nblocks,1,1), 'grid': (nphotons//self.nblocks+1,1)}
- positions_gpu = cuda.mem_alloc(np.dtype(gpuarray.vec.float3).itemsize*nphotons)
- fill_float3(np.int32(nphotons), positions_gpu, position, **gpu_kwargs)
- directions_gpu = cuda.mem_alloc(np.dtype(gpuarray.vec.float3).itemsize*nphotons)
+ positions_gpu = gpuarray.empty(nphotons, dtype=gpuarray.vec.float3)
+ positions_gpu.fill(position)
+ directions_gpu = gpuarray.empty(nphotons, dtype=gpuarray.vec.float3)
fill_uniform_sphere(np.int32(nphotons), rng_states_gpu, directions_gpu, **gpu_kwargs)
- wavelengths_gpu = cuda.mem_alloc(np.dtype(np.float32).itemsize*nphotons)
+ wavelengths_gpu = gpuarray.empty(nphotons, dtype=np.float32)
fill_uniform(np.int32(nphotons), rng_states_gpu, wavelengths_gpu, np.float32(200), np.float32(800), **gpu_kwargs)
- polarizations_gpu = cuda.mem_alloc(np.dtype(gpuarray.vec.float3).itemsize*nphotons)
+ polarizations_gpu = gpuarray.empty(nphotons, dtype=gpuarray.vec.float3)
fill_uniform_sphere(np.int32(nphotons), rng_states_gpu, polarizations_gpu, **gpu_kwargs)
- times_gpu = cuda.mem_alloc(np.dtype(np.float32).itemsize*nphotons)
- fill_float(np.int32(nphotons), times_gpu, np.float32(0), **gpu_kwargs)
- states_gpu = cuda.mem_alloc(np.dtype(np.int32).itemsize*nphotons)
- fill_float(np.int32(nphotons), states_gpu, np.int32(-1), **gpu_kwargs)
- last_hit_triangles_gpu = cuda.mem_alloc(np.dtype(np.int32).itemsize*nphotons)
- fill_float(np.int32(nphotons), last_hit_triangles_gpu, np.int32(-1), **gpu_kwargs)
+ times_gpu = gpuarray.zeros(nphotons, dtype=np.float32)
+ states_gpu = gpuarray.empty(nphotons, dtype=np.int32)
+ states_gpu.fill(-1)
+ last_hit_triangles_gpu = gpuarray.empty(nphotons, dtype=np.int32)
+ last_hit_triangles_gpu.fill(-1)
max_steps = 10
- propagate(np.int32(nphotons), rng_states_gpu, positions_gpu, directions_gpu, wavelengths_gpu, polarizations_gpu, times_gpu, states_gpu, last_hit_triangles_gpu, np.int32(self.geometry.node_map.size-1), np.int32(self.geometry.first_node), np.int32(max_steps), block=(self.nblocks,1,1), grid=(nphotons//self.nblocks+1,1))
+ propagate(np.int32(nphotons), rng_states_gpu, positions_gpu, directions_gpu, wavelengths_gpu, polarizations_gpu, times_gpu, states_gpu, last_hit_triangles_gpu, np.int32(self.geometry.node_map.size-1), np.int32(self.geometry.first_node), np.int32(max_steps), **gpu_kwargs)
reset_earliest_time_int(np.float32(1e9), np.int32(len(earliest_time_int_gpu)), earliest_time_int_gpu, block=(self.nblocks,1,1), grid=(len(earliest_time_int_gpu)//self.nblocks+1,1))
diff --git a/itertoolset.py b/itertoolset.py
index ae4be8a..bd4e45c 100644
--- a/itertoolset.py
+++ b/itertoolset.py
@@ -28,3 +28,11 @@ def unique_everseen(iterable, key=None):
if k not in seen:
seen_add(k)
yield element
+
+def ncycles(iterable, n):
+ "Returns the sequence elements n times"
+ return chain.from_iterable(repeat(tuple(iterable), n))
+
+def take(n, iterable):
+ "Return first n items of the iterable as a list"
+ return list(islice(iterable, n))
diff --git a/make.py b/make.py
index 29f0a04..f30042a 100644
--- a/make.py
+++ b/make.py
@@ -2,7 +2,7 @@ import numpy as np
from geometry import Mesh
from transform import rotate
-def rotate_extrude(x, y, theta=np.pi/32):
+def rotate_extrude(x, y, theta=np.pi/32, remove_duplicate_vertices=True):
x, y, = np.asarray(x), np.asarray(y)
if len(x.shape) != 1 or len(y.shape) != 1 or x.size != y.size:
@@ -34,7 +34,7 @@ def rotate_extrude(x, y, theta=np.pi/32):
triangles[start+step:start+2*step,2] = np.arange(this_slice+1, this_slice+len(points))
triangles[start+step:start+2*step,1] = np.arange(next_slice+1, next_slice+len(points))
- return Mesh(vertices, triangles, remove_duplicate_vertices=True)
+ return Mesh(vertices, triangles, remove_duplicate_vertices=remove_duplicate_vertices)
def cube(size=1):
a = np.sqrt(size**2/2.0)
@@ -49,7 +49,7 @@ def cylinder(radius1=1, radius2=1, height=2, theta=np.pi/32):
def sphere(radius=1, theta=np.pi/32):
profile_angles = np.arange(-np.pi/2, np.pi/2+theta, theta)
- return rotate_extrude(np.cos(profile_angles), np.sin(profile_angles), theta)
+ return rotate_extrude(radius*np.cos(profile_angles), radius*np.sin(profile_angles), theta)
def torus(radius=1, offset=3, phi=np.pi/32, theta=np.pi/32):
profile_angles = np.arange(0, 2*np.pi+phi, phi)
diff --git a/materials/__init__.py b/materials/__init__.py
index d2644ad..f8cf700 100644
--- a/materials/__init__.py
+++ b/materials/__init__.py
@@ -8,25 +8,36 @@ sys.path.append(dir + '/..')
from geometry import Material, Surface
+from pmt_optics import *
+
vacuum = Material('vacuum')
vacuum.set('refractive_index', 1)
vacuum.set('absorption_length', np.finfo(np.float32).max)
vacuum.set('scattering_length', np.finfo(np.float32).max)
lambertian_surface = Surface('lambertian_surface')
-lambertian_surface.set('absorption', 0)
-lambertian_surface.set('reflection_specular', 0)
-lambertian_surface.set('reflection_diffuse', 1)
+lambertian_surface.set('detect', 0)
+lambertian_surface.set('absorb', 0)
+lambertian_surface.set('reflect_specular', 0)
+lambertian_surface.set('reflect_diffuse', 1)
black_surface = Surface('black_surface')
-black_surface.set('absorption', 1)
-black_surface.set('reflection_specular', 0)
-black_surface.set('reflection_diffuse', 0)
+black_surface.set('detect', 0)
+black_surface.set('absorb', 1)
+black_surface.set('reflect_specular', 0)
+black_surface.set('reflect_diffuse', 0)
shiny_surface = Surface('shiny_surface')
-shiny_surface.set('absorption', 0)
-shiny_surface.set('reflection_specular', 1)
-shiny_surface.set('reflection_diffuse', 0)
+shiny_surface.set('detect', 0)
+shiny_surface.set('absorb', 0)
+shiny_surface.set('reflect_specular', 1)
+shiny_surface.set('reflect_diffuse', 0)
+
+glossy_surface = Surface('glossy_surface')
+glossy_surface.set('detect', 0)
+glossy_surface.set('absorb', 0)
+glossy_surface.set('reflect_specular', 0.5)
+glossy_surface.set('reflect_diffuse', 0.5)
db = load(open(dir + '/OPTICS.ratdb'))['OPTICS']
diff --git a/models/portal_gun/Portal_Gun_Barrel_Insert.STL b/models/portal_gun/Portal_Gun_Barrel_Insert.STL
new file mode 100644
index 0000000..277edbc
--- /dev/null
+++ b/models/portal_gun/Portal_Gun_Barrel_Insert.STL
Binary files differ
diff --git a/models/portal_gun/Portal_Gun_Bottom_Cover.STL b/models/portal_gun/Portal_Gun_Bottom_Cover.STL
new file mode 100644
index 0000000..82b54dd
--- /dev/null
+++ b/models/portal_gun/Portal_Gun_Bottom_Cover.STL
Binary files differ
diff --git a/models/portal_gun/Portal_Gun_Center_Tube.STL b/models/portal_gun/Portal_Gun_Center_Tube.STL
new file mode 100644
index 0000000..d35837b
--- /dev/null
+++ b/models/portal_gun/Portal_Gun_Center_Tube.STL
Binary files differ
diff --git a/models/portal_gun/Portal_Gun_Claw.STL b/models/portal_gun/Portal_Gun_Claw.STL
new file mode 100644
index 0000000..310bdae
--- /dev/null
+++ b/models/portal_gun/Portal_Gun_Claw.STL
Binary files differ
diff --git a/models/portal_gun/Portal_Gun_Left_Half.STL b/models/portal_gun/Portal_Gun_Left_Half.STL
new file mode 100644
index 0000000..a96459a
--- /dev/null
+++ b/models/portal_gun/Portal_Gun_Left_Half.STL
Binary files differ
diff --git a/models/portal_gun/Portal_Gun_Right_Half.STL b/models/portal_gun/Portal_Gun_Right_Half.STL
new file mode 100644
index 0000000..3a4c68f
--- /dev/null
+++ b/models/portal_gun/Portal_Gun_Right_Half.STL
Binary files differ
diff --git a/models/portal_gun/Portal_Gun_Top_Claw_Mount.STL b/models/portal_gun/Portal_Gun_Top_Claw_Mount.STL
new file mode 100644
index 0000000..cfe61c7
--- /dev/null
+++ b/models/portal_gun/Portal_Gun_Top_Claw_Mount.STL
Binary files differ
diff --git a/optics.py b/optics.py
new file mode 100644
index 0000000..23a096f
--- /dev/null
+++ b/optics.py
@@ -0,0 +1,148 @@
+import numpy as np
+from geometry import Material, Surface
+
+vacuum = Material('vacuum')
+vacuum.set('refractive_index', 1.0)
+vacuum.set('absorption_length', np.inf)
+vacuum.set('scattering_length', np.inf)
+
+lambertian_surface = Surface('lambertian_surface')
+lambertian_surface.set('reflect_diffuse', 1)
+
+black_surface = Surface('black_surface')
+black_surface.set('absorb', 1)
+
+shiny_surface = Surface('shiny_surface')
+shiny_surface.set('reflect_specular', 1)
+
+glossy_surface = Surface('glossy_surface')
+glossy_surface.set('reflect_diffuse', 0.5)
+glossy_surface.set('reflect_specular', 0.5)
+
+# r7081hqe photocathode material surface
+# source: hamamatsu supplied datasheet for r7081hqe pmt serial number zd0062
+r7081hqe_photocathode = Surface('r7081hqe_photocathode')
+r7081hqe_photocathode.detect = \
+ np.array([(260.0, 0.00),
+ (270.0, 0.04), (280.0, 0.07), (290.0, 0.77), (300.0, 4.57),
+ (310.0, 11.80), (320.0, 17.70), (330.0, 23.50), (340.0, 27.54),
+ (350.0, 30.52), (360.0, 31.60), (370.0, 31.90), (380.0, 32.20),
+ (390.0, 32.00), (400.0, 31.80), (410.0, 30.80), (420.0, 30.16),
+ (430.0, 29.24), (440.0, 28.31), (450.0, 27.41), (460.0, 26.25),
+ (470.0, 24.90), (480.0, 23.05), (490.0, 21.58), (500.0, 19.94),
+ (510.0, 18.48), (520.0, 17.01), (530.0, 15.34), (540.0, 12.93),
+ (550.0, 10.17), (560.0, 7.86), (570.0, 6.23), (580.0, 5.07),
+ (590.0, 4.03), (600.0, 3.18), (610.0, 2.38), (620.0, 1.72),
+ (630.0, 0.95), (640.0, 0.71), (650.0, 0.44), (660.0, 0.25),
+ (670.0, 0.14), (680.0, 0.07), (690.0, 0.03), (700.0, 0.02),
+ (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])
+
+# water data comes from 'lightwater_sno' material in the SNO+ optics database
+water = Material('water')
+water.absorption_length = \
+ np.array([[ 200. , 57.51539993],
+ [ 220. , 64.22219849],
+ [ 240. , 72.6996994 ],
+ [ 260. , 83.75559998],
+ [ 280. , 98.77729797],
+ [ 300. , 120.36499786],
+ [ 320. , 154.0269928 ],
+ [ 340. , 213.82899475],
+ [ 360. , 349.5369873 ],
+ [ 380. , 105.87799835],
+ [ 400. , 50.35989761],
+ [ 420. , 32.56269836],
+ [ 440. , 26.70409966],
+ [ 460. , 22.63209915],
+ [ 480. , 19.63769913],
+ [ 500. , 17.34300041],
+ [ 520. , 11.84370041],
+ [ 540. , 8.99226952],
+ [ 560. , 7.24743032],
+ [ 580. , 6.06968021],
+ [ 600. , 5.22121 ],
+ [ 620. , 4.58085012],
+ [ 640. , 4.08041 ],
+ [ 660. , 3.67853999],
+ [ 680. , 3.3487401 ],
+ [ 700. , 3.07319999],
+ [ 720. , 2.83956003],
+ [ 740. , 2.63893986],
+ [ 760. , 2.46479011],
+ [ 780. , 2.31220984],
+ [ 800. , 2.1774199 ]], dtype=np.float32)
+water.scattering_length = \
+ np.array([[ 200. , 11.36030006],
+ [ 220. , 16.63280106],
+ [ 240. , 23.55719948],
+ [ 260. , 32.44709778],
+ [ 280. , 43.64310074],
+ [ 300. , 57.51350021],
+ [ 320. , 74.45359802],
+ [ 340. , 94.88600159],
+ [ 360. , 119.26100159],
+ [ 380. , 148.05499268],
+ [ 400. , 181.77200317],
+ [ 420. , 220.94500732],
+ [ 440. , 266.13299561],
+ [ 460. , 317.92098999],
+ [ 480. , 376.92300415],
+ [ 500. , 443.78100586],
+ [ 520. , 519.16101074],
+ [ 540. , 603.75897217],
+ [ 560. , 698.29797363],
+ [ 580. , 803.52697754],
+ [ 600. , 920.22399902],
+ [ 620. , 1049.18994141],
+ [ 640. , 1191.27001953],
+ [ 660. , 1347.30004883],
+ [ 680. , 1518.18005371],
+ [ 700. , 1704.82995605],
+ [ 720. , 1908.18005371],
+ [ 740. , 2129.19995117],
+ [ 760. , 2368.87988281],
+ [ 780. , 2628.25 ],
+ [ 800. , 2908.36010742]], dtype=np.float32)
+water.refractive_index = \
+ np.array([[ 200. , 1.41614997],
+ [ 220. , 1.39726996],
+ [ 240. , 1.38395 ],
+ [ 260. , 1.37414002],
+ [ 280. , 1.36667001],
+ [ 300. , 1.36082006],
+ [ 320. , 1.35615003],
+ [ 340. , 1.35232997],
+ [ 360. , 1.34915996],
+ [ 380. , 1.34650004],
+ [ 400. , 1.34423006],
+ [ 420. , 1.34227002],
+ [ 440. , 1.34057999],
+ [ 460. , 1.33908999],
+ [ 480. , 1.33778 ],
+ [ 500. , 1.33660996],
+ [ 520. , 1.33556998],
+ [ 540. , 1.33463001],
+ [ 560. , 1.33378005],
+ [ 580. , 1.33300996],
+ [ 600. , 1.33230996],
+ [ 620. , 1.33167005],
+ [ 640. , 1.33107996],
+ [ 660. , 1.33053994],
+ [ 680. , 1.33003998],
+ [ 700. , 1.32957006],
+ [ 720. , 1.32913995],
+ [ 740. , 1.32874 ],
+ [ 760. , 1.32835996],
+ [ 780. , 1.32800996],
+ [ 800. , 1.32767999]], dtype=np.float32)
+
+# glass data comes from 'glass_sno' material in SNO+ optics database
+glass = Material('glass')
+glass.set('refractive_index', 1.49)
+glass.absorption_length = \
+ np.array([(200, 0.1e-6), (300, 0.1e-6), (330, 1.0), (500, 2.0), (600, 1.0), (770, 0.5), (800, 0.1e-6)])
+glass.set('scattering_length', np.inf)
diff --git a/profiles/hamamatsu_12inch.py b/profiles/hamamatsu_12inch.py
new file mode 100644
index 0000000..c52c029
--- /dev/null
+++ b/profiles/hamamatsu_12inch.py
@@ -0,0 +1,112 @@
+import os
+import sys
+
+dir = os.path.split(os.path.realpath(__file__))[0]
+sys.path.append(dir + '/..')
+
+import numpy as np
+from geometry import Solid
+from make import *
+from view import *
+import matplotlib.pyplot as plt
+from materials import *
+
+f = open('hamamatsu_12inch.txt')
+
+points = []
+for line in f:
+ try:
+ points.append([float(s) for s in line.split(',')])
+ except ValueError:
+ pass
+
+points = np.array(points)
+
+# slice profile in half
+points = points[points[:,0] < 0]
+points[:,0] = -points[:,0]
+# order points from base to face
+points = points[np.argsort(points[:,1])]
+# set x coordinate to 0.0 for first and last points along the profile
+# so that the mesh is closed
+points[0,0] = 0.0
+points[-1,0] = 0.0
+# convert mm -> m
+points /= 1000.0
+
+def offset(points, x):
+ points = np.array([points[0] - (points[1] - points[0])] + list(points) + [points[-1] - (points[-2] - points[-1])])
+
+ offset_points = []
+ for i in range(1,len(points)-1):
+ print '%i' % i
+
+ v1 = np.cross(points[i]-points[i-1], (0,0,1))[:2]
+ v1 /= np.linalg.norm(v1)
+ v1 *= x
+
+ a = points[i-1] + v1
+ b = points[i] + v1
+
+ v2 = np.cross(points[i+1]-points[i], (0,0,1))[:2]
+ v2 /= np.linalg.norm(v2)
+ v2 *= x
+
+ c = points[i] + v2
+ d = points[i+1] + v2
+
+ m = np.empty((2,2))
+ m[:,0] = b-a
+ m[:,1] = c-d
+
+ try:
+ j = np.linalg.solve(m, c-a)[0]
+ except np.linalg.linalg.LinAlgError as e:
+ print e
+ offset_points.append(b)
+ continue
+
+ offset_points.append((a + j*(b-a)))
+
+ return np.array(offset_points)
+
+fig = plt.figure()
+ax = fig.add_subplot(111, aspect='equal')
+
+offset_points = offset(points, -.003)
+
+plt.plot(points[:,0], points[:,1], 'b-', offset_points[:,0], offset_points[:,1], 'g-')
+plt.show()
+
+pmt_outer_mesh = rotate_extrude(points[:,0], points[:,1], np.pi/8)
+pmt_inner_mesh = rotate_extrude(offset_points[:,0], offset_points[:,1], np.pi/8)
+
+# cutaway view
+pmt_outer_mesh.triangles = pmt_outer_mesh.triangles[np.mean(pmt_outer_mesh[:], axis=1)[:,0] > 0]
+
+#pmt_outer_mesh.triangles = pmt_outer_mesh.triangles[np.mean(pmt_outer_mesh[:], axis=1)[:,1] > -1e-3]
+#pmt_inner_mesh.triangles = pmt_inner_mesh.triangles[np.mean(pmt_inner_mesh[:], axis=1)[:,1] > -1e-3]
+
+pmt_outer_solid = Solid(pmt_outer_mesh)
+
+
+pmt_inner_solid = Solid(pmt_inner_mesh, color=0x00ff00)
+
+def get_lc_profile(radii, a, b, d, rmin, rmax):
+ c = -b*np.sqrt(1 - (rmin-d)**2/a**2)
+
+ return -c - b*np.sqrt(1-(radii-d)**2/a**2)
+
+lc_radii = np.linspace(152.4e-3, 209.672e-3 + 152.4e-3 - 127e-3, 10)
+lc_profile = get_lc_profile(lc_radii, 165.97e-3, 584.525e-3, 95.48e-3, 152.4e-3, 209.672e-3 + 152.4e-3 - 127e-3)
+
+face_points = points[points[:,1] > -1e-3]
+
+lc_offset = np.interp(lc_radii[0], list(reversed(face_points[:,0])), list(reversed(face_points[:,1])))
+
+lc_mesh = rotate_extrude(lc_radii, lc_profile + lc_offset, np.pi/8)
+lc_solid = Solid(lc_mesh, color=0xff0000, surface=shiny_surface)
+
+pmt_solid_lc = pmt_inner_solid + pmt_outer_solid + lc_solid
+
+view(pmt_solid_lc)
diff --git a/profiles/hamamatsu_12inch.txt b/profiles/hamamatsu_12inch.txt
new file mode 100644
index 0000000..e896b94
--- /dev/null
+++ b/profiles/hamamatsu_12inch.txt
@@ -0,0 +1,55 @@
+#DataThief /Users/stan/Downloads/./Dimension_of_12inch_bulb_100826.png Monday 16-May-2011 12:22:34 PM
+109.1514, 83.9786
+123.0527, 70.9583
+134.8946, 56.9466
+145.1919, 40.944
+150.8554, 23.4594
+152.4, 3.9922
+150.8554, -15.4649
+145.7067, -32.914
+136.4392, -49.3517
+123.0527, -64.2787
+108.6364, -75.2102
+92.1608, -86.1349
+76.2, -97.0612
+66.4176, -112.0001
+63.3284, -130.9531
+63.3284, -150.9143
+63.3284, -170.3764
+56.6351, -187.8206
+44.7932, -202.7527
+28.8324, -211.1839
+10.8122, -212.6219
+-9.2676, -213.5543
+-28.8324, -210.995
+-45.823, -201.9569
+-57.6649, -186.9473
+-63.8432, -169.96
+-63.8432, -150.4979
+-63.8432, -130.5367
+-66.4176, -111.0661
+-76.7149, -96.0615
+-93.1905, -85.0289
+-108.6364, -74.4987
+-123.0527, -63.4727
+-135.4095, -49.9586
+-144.677, -33.9593
+-150.3405, -16.4747
+-152.9149, 2.4968
+-151.3703, 21.4549
+-146.2216, 39.4031
+-136.4392, 55.34
+-124.0824, 70.2704
+-109.6662, 83.198
+-93.7054, 94.1244
+-76.7149, 103.0513
+-59.2095, 109.4813
+-40.6743, 113.9119
+-21.1095, 117.3411
+-1.0297, 118.2734
+19.05, 117.2096
+38.6149, 114.1513
+57.15, 109.5994
+75.1703, 103.053
+92.1608, 94.5138
+104.0027, 88.4867
diff --git a/profiles/test.py b/profiles/test.py
new file mode 100644
index 0000000..2ff3168
--- /dev/null
+++ b/profiles/test.py
@@ -0,0 +1,3 @@
+hello world
+
+hello world
diff --git a/render.py b/render.py
new file mode 100755
index 0000000..2a4e0d5
--- /dev/null
+++ b/render.py
@@ -0,0 +1,468 @@
+#!/usr/bin/env python
+import os
+import sys
+import numpy as np
+import pickle
+import tempfile
+
+import pygame
+from pygame.locals import *
+
+import src
+from camera import *
+from geometry import *
+from transform import *
+from optics import *
+
+import time
+
+from pycuda import autoinit
+from pycuda.compiler import SourceModule
+from pycuda import gpuarray
+from pycuda.characterize import sizeof
+
+def to_geometry(viewable, bits):
+ if isinstance(viewable, Geometry):
+ geometry = viewable
+ geometry.build(bits)
+ elif isinstance(viewable, Solid):
+ geometry = Geometry()
+ geometry.add_solid(viewable)
+ geometry.build(bits)
+ elif isinstance(viewable, Mesh):
+ geometry = Geometry()
+ geometry.add_solid(Solid(viewable, vacuum, vacuum, surface=lambertian_surface))
+ geometry.build(bits)
+ else:
+ raise Exception("can't convert %s to a geometry" % viewable)
+
+ return geometry
+
+def box(lower_bound, upper_bound):
+ dx, dy, dz = upper_bound - lower_bound
+
+ vertices = np.array([lower_bound,
+ lower_bound + (dx,0,0),
+ lower_bound + (dx,dy,0),
+ lower_bound + (0,dy,0),
+ lower_bound + (0,0,dz),
+ lower_bound + (dx,0,dz),
+ lower_bound + (dx,dy,dz),
+ lower_bound + (0,dy,dz)])
+
+ triangles = np.empty((12,3), dtype=np.int32)
+ # bottom
+ triangles[0] = (1,3,2)
+ triangles[1] = (1,4,3)
+ # top
+ triangles[2] = (5,7,6)
+ triangles[3] = (5,8,7)
+ # left
+ triangles[4] = (5,1,2)
+ triangles[5] = (5,2,6)
+ # right
+ triangles[6] = (3,4,8)
+ triangles[7] = (3,8,7)
+ # front
+ triangles[8] = (2,3,7)
+ triangles[9] = (2,7,6)
+ # back
+ triangles[10] = (1,5,8)
+ triangles[11] = (1,8,4)
+
+ triangles -= 1
+
+ return Mesh(vertices, triangles)
+
+def bvh_mesh(geometry, layer):
+ lower_bounds = geometry.lower_bounds[geometry.layers == layer]
+ upper_bounds = geometry.upper_bounds[geometry.layers == layer]
+
+ if len(lower_bounds) == 0 or len(upper_bounds) == 0:
+ raise Exception('no nodes at layer %i' % layer)
+
+ mesh = box(lower_bounds[0], upper_bounds[0])
+
+ for lower_bound, upper_bound in zip(lower_bounds, upper_bounds)[1:]:
+ mesh += box(lower_bound, upper_bound)
+
+ return mesh
+
+def view(viewable, name='', bits=8, load_bvh=False, camera_view=None, make_movie=False):
+ """
+ Render `viewable` in a pygame window.
+
+ Movement:
+ - zoom: scroll the mouse wheel
+ - rotate: click and drag the mouse
+ - move: shift+click and drag the mouse
+ """
+
+ geometry = to_geometry(viewable, bits)
+
+ if load_bvh:
+ print 'loading bvh...'
+ bvhg = []
+
+ bvhg.append(geometry)
+
+ for layer in sorted(np.unique(geometry.layers)):
+ print 'building layer %i' % layer
+ bvhg.append(to_geometry(bvh_mesh(geometry, layer), bits))
+
+ lower_bound = np.array([np.min(geometry.mesh[:][:,:,0]),
+ np.min(geometry.mesh[:][:,:,1]),
+ np.min(geometry.mesh[:][:,:,2])])
+
+ upper_bound = np.array([np.max(geometry.mesh[:][:,:,0]),
+ np.max(geometry.mesh[:][:,:,1]),
+ np.max(geometry.mesh[:][:,:,2])])
+
+ scale = np.linalg.norm(upper_bound-lower_bound)
+
+ print 'device %s' % autoinit.device.name()
+
+ module = SourceModule(src.kernel, options=['-I' + src.dir],
+ no_extern_c=True, cache_dir=False)
+ #geometry.load(module, color=True)
+ geometry.load(module, color=False)
+ #cuda_raytrace = module.get_function('ray_trace')
+ cuda_build_rgb_lookup = module.get_function('build_rgb_lookup')
+ cuda_render = module.get_function('render')
+ cuda_rotate = module.get_function('rotate')
+ cuda_translate = module.get_function('translate')
+
+ pygame.init()
+ size = width, height = 800, 600
+
+ init_rng = module.get_function('init_rng')
+ rng_states_gpu = cuda.mem_alloc(len(geometry.mesh.triangles)*sizeof('curandStateXORWOW', '#include <curand_kernel.h>'))
+ init_rng(np.int32(len(geometry.mesh.triangles)), rng_states_gpu, np.int32(0), np.int32(0), block=(64,1,1), grid=(len(geometry.mesh.triangles)//64+1,1))
+
+ diagonal = np.linalg.norm(upper_bound-lower_bound)
+ #source_position = [0, diagonal*1.75, (lower_bound[2]+upper_bound[2])/2]
+ source_position = [0, 0, upper_bound[2]+5.0]
+
+ source_positions_gpu = gpuarray.empty(len(geometry.mesh.triangles), dtype=gpuarray.vec.float3)
+ source_positions_gpu.fill(gpuarray.vec.make_float3(*source_position))
+
+ source_directions = np.mean(geometry.mesh[:], axis=1) - source_position
+ source_directions_float3 = np.empty(source_positions_gpu.size, dtype=gpuarray.vec.float3)
+ source_directions_float3['x'] = source_directions[:,0]
+ source_directions_float3['y'] = source_directions[:,1]
+ source_directions_float3['z'] = source_directions[:,2]
+ source_directions_gpu = cuda.to_device(source_directions_float3)
+
+ if (len(geometry.mesh.triangles) < source_positions_gpu.size):
+ print width*height
+ print source_positions_gpu.size
+ raise Exception('not enough rng_states')
+
+ rgb_lookup_gpu = gpuarray.zeros(source_positions_gpu.size, dtype=gpuarray.vec.float3)
+
+ max_steps = 10
+ rgb_runs = 1000
+
+ cuda_build_rgb_lookup(np.int32(source_positions_gpu.size), rng_states_gpu, source_positions_gpu, source_directions_gpu, np.int32(geometry.node_map.size-1), np.int32(geometry.first_node), rgb_lookup_gpu, np.int32(rgb_runs), np.int32(max_steps), block=(64,1,1), grid=(source_positions_gpu.size//64+1,1))
+
+ rgb_lookup = rgb_lookup_gpu.get()
+ rgb_lookup['x'] /= rgb_runs
+ rgb_lookup['y'] /= rgb_runs
+ rgb_lookup['z'] /= rgb_runs
+ rgb_lookup_gpu = cuda.to_device(rgb_lookup)
+
+
+ print rgb_lookup
+
+
+
+
+
+
+
+
+
+
+
+
+ if camera_view is None:
+ camera = Camera(size)
+
+ diagonal = np.linalg.norm(upper_bound-lower_bound)
+
+ 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)
+
+ origins, directions = camera.get_rays()
+
+ origins_float3 = np.empty(origins.shape[0], dtype=gpuarray.vec.float3)
+ origins_float3['x'] = origins[:,0]
+ origins_float3['y'] = origins[:,1]
+ origins_float3['z'] = origins[:,2]
+
+ directions_float3 = np.empty(directions.shape[0], dtype=gpuarray.vec.float3)
+ directions_float3['x'] = directions[:,0]
+ directions_float3['y'] = directions[:,1]
+ directions_float3['z'] = directions[:,2]
+ else:
+ f = open(camera_view, 'rb')
+ origins_float3 = pickle.load(f)
+ directions_float3 = pickle.load(f)
+ point = pickle.load(f)
+ axis1 = pickle.load(f)
+ axis2 = pickle.load(f)
+ f.close()
+
+ origins_gpu = cuda.to_device(origins_float3)
+ directions_gpu = cuda.to_device(directions_float3)
+
+ pixels = np.empty(width*height, dtype=np.int32)
+ pixels_gpu = cuda.to_device(pixels)
+
+ nruns = 5
+
+
+
+ screen = pygame.display.set_mode(size)
+ pygame.display.set_caption(name)
+
+
+
+
+
+
+ nblocks = 64
+
+ gpu_kwargs = {'block': (nblocks,1,1), 'grid':(pixels.size/nblocks+1,1)}
+
+ if make_movie:
+ tempdir = tempfile.mkdtemp()
+
+ def screenshot(dir=''):
+ if name == '':
+ root, ext = 'screenshot', 'png'
+ else:
+ root, ext = name, 'png'
+
+ filename = os.path.join(dir,'.'.join([root, ext]))
+
+ i = 1
+ while os.path.exists(filename):
+ filename = os.path.join(dir,'.'.join([root + str(i), ext]))
+ i += 1
+
+ pygame.image.save(screen, filename)
+ print 'image saved to %s' % filename
+
+ def render():
+ """Render the mesh and display to screen."""
+ t0 = time.time()
+ cuda_render(np.int32(pixels.size), rng_states_gpu, origins_gpu, directions_gpu, np.int32(geometry.node_map.size-1), np.int32(geometry.first_node), rgb_lookup_gpu, np.int32(nruns), pixels_gpu, np.int32(max_steps), **gpu_kwargs)
+ cuda.Context.synchronize()
+ elapsed = time.time() - t0
+
+ print 'elapsed %f sec' % elapsed
+
+ cuda.memcpy_dtoh(pixels, pixels_gpu)
+ pygame.surfarray.blit_array(screen, pixels.reshape(size))
+ pygame.display.flip()
+
+ if make_movie:
+ screenshot(tempdir)
+
+ render()
+
+ 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 = scale*np.cross(axis1,axis2)/10.0
+
+ cuda_translate(np.int32(pixels.size), origins_gpu, gpuarray.vec.make_float3(*v), **gpu_kwargs)
+
+ point += v
+
+ render()
+
+ if event.button == 5:
+ v = -scale*np.cross(axis1,axis2)/10.0
+
+ cuda_translate(np.int32(pixels.size), origins_gpu, gpuarray.vec.make_float3(*v), **gpu_kwargs)
+
+ point += v
+
+ render()
+
+ 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]*axis1 + movement[1]*axis2
+ mouse_direction /= np.linalg.norm(mouse_direction)
+
+ if shift:
+ v = mouse_direction*scale*length/float(width)
+
+ cuda_translate(np.int32(pixels.size), origins_gpu, gpuarray.vec.make_float3(*v), **gpu_kwargs)
+
+ point += v
+
+ render()
+ else:
+ phi = np.float32(2*np.pi*length/float(width))
+ n = rotate(mouse_direction, np.pi/2, \
+ -np.cross(axis1,axis2))
+
+ cuda_rotate(np.int32(pixels.size), origins_gpu, phi, gpuarray.vec.make_float3(*n), **gpu_kwargs)
+
+ cuda_rotate(np.int32(pixels.size), directions_gpu, phi, gpuarray.vec.make_float3(*n), **gpu_kwargs)
+
+ point = rotate(point, phi, n)
+ axis1 = rotate(axis1, phi, n)
+ axis2 = rotate(axis2, phi, n)
+
+ render()
+
+ if event.type == KEYDOWN:
+ 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:
+ try:
+ if current_layer+1 >= len(bvhg):
+ raise IndexError
+
+ geometry = bvhg[current_layer+1]
+ current_layer += 1
+
+ geometry.load(module, color=True)
+ render()
+ except IndexError:
+ print 'no further layers to view'
+
+ if event.key == K_PAGEDOWN:
+ try:
+ if current_layer-1 < 0:
+ raise IndexError
+
+ geometry = bvhg[current_layer-1]
+ current_layer -= 1
+
+ geometry.load(module, color=True)
+ render()
+ except IndexError:
+ print 'no further layers to view'
+
+ if event.key == K_F11:
+ if name == '':
+ root, ext = 'camera_view', 'pkl'
+ else:
+ root, ext = name + '_camera_view', 'pkl'
+
+ filename = '.'.join([root, ext])
+
+ i = 1
+ while os.path.exists(filename):
+ filename = '.'.join([root + str(i), ext])
+ i += 1
+
+ f = open(filename, 'wb')
+ cuda.memcpy_dtoh(origins_float3, origins_gpu)
+ cuda.memcpy_dtoh(directions_float3, directions_gpu)
+ pickle.dump(origins_float3, f)
+ pickle.dump(directions_float3, f)
+ pickle.dump(point, f)
+ pickle.dump(axis1, f)
+ pickle.dump(axis2, f)
+ f.close()
+
+ print 'camera view saved to %s' % filename
+
+ if event.key == K_F12:
+ screenshot()
+
+ if event.type == KEYUP:
+ if event.key == K_LSHIFT or event.key == K_RSHIFT:
+ shift = False
+
+ pygame.display.quit()
+
+ if make_movie:
+ if name == '':
+ root, ext = 'movie', 'tgz'
+ else:
+ root, ext = name + '_movie', 'tgz'
+
+ filename = '.'.join([root, ext])
+
+ i=1
+ while os.path.exists(filename):
+ filename = '.'.join([root + str(i), ext])
+ i += 1
+
+ from subprocess import call
+ call(['tar', '-czvf', filename, tempdir])
+
+ os.removedirs(tempdir)
+
+if __name__ == '__main__':
+ import optparse
+
+ from stl import mesh_from_stl
+
+ import solids
+ import detectors
+
+ 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('-v', '--camera_view', dest='camera_view',
+ help='load a camera view', default=None)
+ parser.add_option('-m', '--movie', action='store_true', dest='movie',
+ help='create a movie', default=False)
+ options, args = parser.parse_args()
+
+ if len(args) < 1:
+ sys.exit(parser.format_help())
+
+ head, tail = os.path.split(args[0])
+ root, ext = os.path.splitext(tail)
+
+ if ext.lower() == '.stl':
+ view(mesh_from_stl(args[0]), root, options.bits, options.load_bvh, options.camera_view, options.movie)
+ else:
+ import inspect
+
+ members = dict(inspect.getmembers(detectors) + inspect.getmembers(solids))
+
+ if args[0] in members:
+ view(members[args[0]], args[0], options.bits, options.load_bvh, options.camera_view, options.movie)
+ else:
+ sys.exit("couldn't find object %s" % args[0])
diff --git a/scenes/checkerboard.py b/scenes/checkerboard.py
new file mode 100644
index 0000000..2852348
--- /dev/null
+++ b/scenes/checkerboard.py
@@ -0,0 +1,58 @@
+import os
+import sys
+
+dir = os.path.split(os.path.realpath(__file__))[0]
+sys.path.append(dir + '/..')
+
+import numpy as np
+from itertools import product
+from itertoolset import ncycles
+from geometry import *
+from optics import *
+from make import *
+
+checkers_per_side = 10
+squares_per_checker = 10
+
+
+x = np.linspace(-5.0, 5.0, checkers_per_side*squares_per_checker+1)
+y = np.linspace(-5.0, 5.0, checkers_per_side*squares_per_checker+1)
+
+vertices = np.array(tuple(product(x,y,[0])))
+
+triangles = []
+for j in range(y.size-1):
+ for i in range(x.size-1):
+ triangles.append([j*len(x)+i, (j+1)*len(x)+i,(j+1)*len(x)+i+1])
+ triangles.append([j*len(x)+i, j*len(x)+i+1,(j+1)*len(x)+i+1])
+
+checkerboard = Mesh(vertices, triangles, remove_duplicate_vertices=True)
+checker_color_line1 = take(checkers_per_side*squares_per_checker*2, cycle([0]*2*squares_per_checker + [0xffffff]*2*squares_per_checker))*squares_per_checker
+checker_color_line2 = take(checkers_per_side*squares_per_checker*2, cycle([0xffffff]*2*squares_per_checker + [0]*2*squares_per_checker))*squares_per_checker
+checker_color = take(len(checkerboard.triangles), cycle(checker_color_line1 + checker_color_line2))
+
+print len(checkerboard.triangles)
+print len(checker_color)
+
+checker_surface_line1 = take(checkers_per_side*squares_per_checker*2, cycle([black_surface]*2*squares_per_checker + [lambertian_surface]*2*squares_per_checker))*squares_per_checker
+checker_surface_line2 = take(checkers_per_side*squares_per_checker*2, cycle([lambertian_surface]*2*squares_per_checker + [black_surface]*2*squares_per_checker))*squares_per_checker
+checker_surface = take(len(checkerboard.triangles), cycle(checker_surface_line1 + checker_surface_line2))
+
+checkerboard_solid = Solid(checkerboard, vacuum, vacuum, surface=checker_surface, color=checker_color)
+
+print len(checker_surface)
+
+water_sphere_solid = Solid(sphere(theta=np.pi/64), water, vacuum)
+mirror_sphere_solid = Solid(sphere(theta=np.pi/64), vacuum, vacuum, surface=shiny_surface)
+diffuse_sphere_solid = Solid(sphere(theta=np.pi/64), vacuum, vacuum, surface=lambertian_surface)
+
+checkerboard_scene = Geometry()
+checkerboard_scene.add_solid(checkerboard_solid, displacement=(0,0,-1.5))
+checkerboard_scene.add_solid(water_sphere_solid, displacement=(2.0,-2.0,0))
+checkerboard_scene.add_solid(mirror_sphere_solid, displacement=(-2.0,-2.0,0))
+checkerboard_scene.add_solid(diffuse_sphere_solid, displacement=(0.0,2.0,0))
+
+#from view import *
+from render import *
+
+view(checkerboard_scene)
diff --git a/solids/r11708.py b/solids/r11708.py
index a0476d0..d1a5514 100644
--- a/solids/r11708.py
+++ b/solids/r11708.py
@@ -20,7 +20,7 @@ inner_color[photocathode_triangles] = 0xff0000
inner_color[~photocathode_triangles] = 0x00ff00
inner_surface = np.empty(len(r11708_inner_mesh.triangles), np.object)
-inner_surface[photocathode_triangles] = black_surface
+inner_surface[photocathode_triangles] = r7081hqe_photocathode
inner_surface[~photocathode_triangles] = shiny_surface
r11708_inner_solid = Solid(r11708_inner_mesh, vacuum, glass, inner_surface, color=inner_color)
diff --git a/solids/r11708_cut.py b/solids/r11708_cut.py
index a9ffeba..14cfff0 100644
--- a/solids/r11708_cut.py
+++ b/solids/r11708_cut.py
@@ -20,7 +20,7 @@ inner_color[photocathode_triangles] = 0xff0000
inner_color[~photocathode_triangles] = 0x00ff00
inner_surface = np.empty(len(r11708_inner_mesh.triangles), np.object)
-inner_surface[photocathode_triangles] = black_surface
+inner_surface[photocathode_triangles] = r7081hqe_photocathode
inner_surface[~photocathode_triangles] = shiny_surface
r11708_inner_solid = Solid(r11708_inner_mesh, vacuum, glass, inner_surface, color=inner_color)
diff --git a/src/kernel.cu b/src/kernel.cu
index f257c6c..aae6e95 100644
--- a/src/kernel.cu
+++ b/src/kernel.cu
@@ -14,10 +14,12 @@
enum
{
+ DIFFUSE_HIT = -3,
DEBUG = -2,
INIT = -1,
NO_HIT,
BULK_ABSORB,
+ SURFACE_DETECT,
SURFACE_ABSORB
};
@@ -144,28 +146,235 @@ __device__ int intersect_mesh(const float3 &origin, const float3& direction, con
return triangle_index;
}
-extern "C"
+__device__ void myAtomicAdd(float *addr, float data)
{
+ while (data)
+ data = atomicExch(addr, data+atomicExch(addr, 0.0f));
+}
+
+
-__global__ void fill_float(int nthreads, float *a, float value)
-{
- int id = blockIdx.x*blockDim.x + threadIdx.x;
- if (id >= nthreads)
- return;
- a[id] = value;
-}
-__global__ void fill_float3(int nthreads, float3 *a, float3 value)
+
+__device__ int to_diffuse(curandState &rng, float3 position, float3 direction, float wavelength, float3 polarization, int start_node, int first_node, int max_steps)
{
- int id = blockIdx.x*blockDim.x + threadIdx.x;
+ int last_hit_triangle = -1;
+
+ int steps = 0;
+ while (steps < max_steps)
+ {
+ steps++;
+
+ float distance;
+
+ last_hit_triangle = intersect_mesh(position, direction, start_node, first_node, distance, last_hit_triangle);
+
+ if (last_hit_triangle == -1)
+ return last_hit_triangle;
+
+ uint4 triangle_data = triangles[last_hit_triangle];
+
+ float3 v0 = vertices[triangle_data.x];
+ float3 v1 = vertices[triangle_data.y];
+ float3 v2 = vertices[triangle_data.z];
+
+ int material_in_index = convert(0xFF & (triangle_data.w >> 24));
+ int material_out_index = convert(0xFF & (triangle_data.w >> 16));
+ int surface_index = convert(0xFF & (triangle_data.w >> 8));
+
+ float3 surface_normal = cross(v1-v0,v2-v1);
+ surface_normal /= norm(surface_normal);
+
+ Material material1, material2;
+ if (dot(surface_normal,-direction) > 0.0f)
+ {
+ // outside to inside
+ material1 = materials[material_out_index];
+ material2 = materials[material_in_index];
+ }
+ else
+ {
+ // inside to outside
+ material1 = materials[material_in_index];
+ material2 = materials[material_out_index];
+ surface_normal = -surface_normal;
+ }
+
+ float refractive_index1 = interp_property(wavelength, material1.refractive_index);
+ float refractive_index2 = interp_property(wavelength, material2.refractive_index);
+
+ float absorption_length = interp_property(wavelength, material1.absorption_length);
+ float scattering_length = interp_property(wavelength, material1.scattering_length);
+
+ float absorption_distance = -absorption_length*logf(curand_uniform(&rng));
+ float scattering_distance = -scattering_length*logf(curand_uniform(&rng));
+
+ if (absorption_distance <= scattering_distance)
+ {
+ if (absorption_distance <= distance)
+ {
+ return -1;
+ } // photon is absorbed in material1
+ }
+ else
+ {
+ if (scattering_distance <= distance)
+ {
+ //time += scattering_distance/(SPEED_OF_LIGHT/refractive_index1);
+ position += scattering_distance*direction;
+
+ float theta,y;
+ while (true)
+ {
+ y = curand_uniform(&rng);
+ theta = uniform(&rng, 0, 2*PI);
+
+ if (y < powf(cosf(theta),2))
+ break;
+ }
+
+ float phi = uniform(&rng, 0, 2*PI);
+
+ float3 b = cross(polarization, direction);
+ float3 c = polarization;
+
+ direction = rotate(direction, theta, b);
+ direction = rotate(direction, phi, c);
+ polarization = rotate(polarization, theta, b);
+ polarization = rotate(polarization, phi, c);
+
+ last_hit_triangle = -1;
+
+ continue;
+ } // photon is scattered in material1
+
+ } // if scattering_distance < absorption_distance
+
+ position += distance*direction;
+ //time += distance/(SPEED_OF_LIGHT/refractive_index1);
+
+ // p is normal to the plane of incidence
+ float3 p = cross(direction, surface_normal);
+ p /= norm(p);
+
+ float normal_coefficient = dot(polarization, p);
+ float normal_probability = normal_coefficient*normal_coefficient;
+
+ float incident_angle = acosf(dot(surface_normal,-direction));
+
+ if (surface_index != -1)
+ {
+ Surface surface = surfaces[surface_index];
+
+ float detect = interp_property(wavelength, surface.detect);
+ float absorb = interp_property(wavelength, surface.absorb);
+ float reflect_diffuse = interp_property(wavelength, surface.reflect_diffuse);
+ float reflect_specular = interp_property(wavelength, surface.reflect_specular);
+
+ // since the surface properties are interpolated linearly,
+ // we are guaranteed that they still sum to one.
+
+ float uniform_sample = curand_uniform(&rng);
+
+ if (uniform_sample < absorb)
+ {
+ // absorb
+ //state = SURFACE_ABSORB;
+ //break;
+ return -1;
+ }
+ else if (uniform_sample < absorb + detect)
+ {
+ // detect
+ //state = SURFACE_DETECT;
+ //break;
+ return -1;
+ }
+ else if (uniform_sample < absorb + detect + reflect_diffuse)
+ {
+ //state = DIFFUSE_HIT;
+ //break;
+ return last_hit_triangle;
+ }
+ else
+ {
+ // specularly reflect
+ direction = rotate(surface_normal, incident_angle, p);
+
+ // polarization ?
+
+ continue;
+ }
+
+ } // surface
+
+ float refracted_angle = asinf(sinf(incident_angle)*refractive_index1/refractive_index2);
+
+ if (curand_uniform(&rng) < normal_probability)
+ {
+
+ float reflection_coefficient = -sinf(incident_angle-refracted_angle)/sinf(incident_angle+refracted_angle);
+ float reflection_probability = reflection_coefficient*reflection_coefficient;
+
+ if ((curand_uniform(&rng) < reflection_probability) || isnan(refracted_angle))
+ {
+ direction = rotate(surface_normal, incident_angle, p);
+ }
+ else
+ {
+ direction = rotate(surface_normal, PI-refracted_angle, p);
+ }
+
+ polarization = p;
+
+ continue;
+ } // photon polarization normal to plane of incidence
+ else
+ {
+ float reflection_coefficient = tanf(incident_angle-refracted_angle)/tanf(incident_angle+refracted_angle);
+ float reflection_probability = reflection_coefficient*reflection_coefficient;
+
+ if ((curand_uniform(&rng) < reflection_probability) || isnan(refracted_angle))
+ {
+ direction = rotate(surface_normal, incident_angle, p);
+ }
+ else
+ {
+ direction = rotate(surface_normal, PI-refracted_angle, p);
+ }
+
+ polarization = cross(p, direction);
+ polarization /= norm(polarization);
+
+ continue;
+ } // photon polarization parallel to surface
+
+ } // while(nsteps < max_steps)
+
+ return -1;
+
+} // to_diffuse
+
+
- if (id >= nthreads)
- return;
- a[id] = value;
-}
+
+
+
+
+
+
+
+
+
+
+
+
+
+extern "C"
+{
__global__ void set_pointer(uint4 *triangle_ptr, float3 *vertex_ptr)
{
@@ -196,6 +405,197 @@ __global__ void rotate(int nthreads, float3 *points, float phi, float3 axis)
points[id] = rotate(points[id], phi, axis);
}
+#define RED_WAVELENGTH 685
+#define BLUE_WAVELENGTH 465
+#define GREEN_WAVELENGTH 545
+
+
+__global__ void build_rgb_lookup(int nthreads, curandState *rng_states, float3 *positions, float3 *directions, int start_node, int first_node, float3 *rgb_lookup, int runs, int max_steps)
+{
+ int id = blockIdx.x*blockDim.x + threadIdx.x;
+
+ if (id >= nthreads)
+ return;
+
+ curandState rng = rng_states[id];
+ float3 position = positions[id];
+ float3 direction = directions[id];
+ direction /= norm(direction);
+ float3 polarization = uniform_sphere(&rng);
+
+ float distance;
+
+ int hit_triangle = intersect_mesh(position, direction, start_node, first_node, distance);
+
+ if (hit_triangle != id)
+ return;
+
+ // note triangles from built geometry mesh are in order
+
+ uint4 triangle_data = triangles[hit_triangle];
+
+ float3 v0 = vertices[triangle_data.x];
+ float3 v1 = vertices[triangle_data.y];
+ float3 v2 = vertices[triangle_data.z];
+
+ float cos_theta = dot(normalize(cross(v1-v0, v2-v1)), -direction);
+
+ if (cos_theta < 0.0f)
+ cos_theta = dot(-normalize(cross(v1-v0, v2-v1)), -direction);
+
+ for (int i=0; i < runs; i++)
+ {
+ hit_triangle = to_diffuse(rng, position, direction, RED_WAVELENGTH, polarization, start_node, first_node, max_steps);
+
+ if (hit_triangle != -1)
+ myAtomicAdd(&rgb_lookup[hit_triangle].x, cos_theta);
+
+ hit_triangle = to_diffuse(rng, position, direction, BLUE_WAVELENGTH, polarization, start_node, first_node, max_steps);
+
+ if (hit_triangle != -1)
+ myAtomicAdd(&rgb_lookup[hit_triangle].y, cos_theta);
+
+ hit_triangle = to_diffuse(rng, position, direction, GREEN_WAVELENGTH, polarization, start_node, first_node, max_steps);
+
+ if (hit_triangle != -1)
+ myAtomicAdd(&rgb_lookup[hit_triangle].z, cos_theta);
+ }
+
+} // build_rgb_lookup
+
+__global__ void render(int nthreads, curandState *rng_states, float3 *positions, float3 *directions, int start_node, int first_node, float3 *rgb_lookup, int runs, int *pixels, int max_steps)
+{
+ int id = blockIdx.x*blockDim.x + threadIdx.x;
+
+ if (id >= nthreads)
+ return;
+
+ curandState rng = rng_states[id];
+ float3 position = positions[id];
+ float3 direction = directions[id];
+ direction /= norm(direction);
+ float3 polarization = uniform_sphere(&rng);
+
+ float3 rgb = make_float3(0.0, 0.0, 0.0);
+
+ int hit_triangle;
+
+ for (int i=0; i < runs; i++)
+ {
+ hit_triangle = to_diffuse(rng, position, direction, RED_WAVELENGTH, polarization, start_node, first_node, max_steps);
+
+ if (hit_triangle != -1)
+ rgb.x += rgb_lookup[hit_triangle].x;
+
+ hit_triangle = to_diffuse(rng, position, direction, BLUE_WAVELENGTH, polarization, start_node, first_node, max_steps);
+
+ if (hit_triangle != -1)
+ rgb.y += rgb_lookup[hit_triangle].y;
+
+ hit_triangle = to_diffuse(rng, position, direction, GREEN_WAVELENGTH, polarization, start_node, first_node, max_steps);
+
+ if (hit_triangle != -1)
+ rgb.z += rgb_lookup[hit_triangle].z;
+ }
+
+ rgb /= runs;
+
+ unsigned int r = floorf(rgb.x*255);
+ unsigned int g = floorf(rgb.y*255);
+ unsigned int b = floorf(rgb.z*255);
+
+ pixels[id] = r << 16 | g << 8 | b;
+
+} // render
+
+#if 0
+
+__global__ void get_triangle(int nthreads, float3 *positions, float3 *directions, int start_node, int first_node, int *triangles)
+{
+ int id = blockIdx.x*blockDim.x + threadIdx.x;
+
+ if (id >= nthreads)
+ return;
+
+ float3 position = positions[id];
+ float3 direction = directions[id];
+ direction /= norm(direction);
+
+ float distance;
+
+ triangles[id] = intersect_mesh(position, direction, start_node, first_node, distance);
+}
+
+__global__ void get_cos_theta(int nthreads, float3 *positions, float3 *directions, int start_node, int first_node, int *triangle, float *cos_thetas)
+{
+ int id = blockIdx.x*blockDim.x + threadIdx.x;
+
+ if (id >= nthreads)
+ return;
+
+ float3 position = positions[id];
+ float3 direction = directions[id];
+ direction /= norm(direction);
+
+ uint4 triangle_data = triangles[triangle[id]];
+
+ float3 v0 = vertices[triangle_data.x];
+ float3 v1 = vertices[triangle_data.y];
+ float3 v2 = vertices[triangle_data.z];
+
+ float cos_theta = dot(normalize(cross(v1-v0, v2-v1)), -direction);
+
+ if (cos_theta < 0.0f)
+ cos_theta = dot(-normalize(cross(v1-v0, v2-v1)), -direction);
+
+ cos_thetas[id] = cos_theta;
+}
+
+#endif
+
+#if 0
+
+__global__ void to_diffuse(int nthreads, curandState *rng_states, float3 *positions, float3 *directions, float *wavelengths, float3 *polarizations, int *last_hit_triangles, int start_node, int first_node, int max_steps)
+{
+ int id = blockIdx.x*blockDim.x + threadIdx.x;
+
+ if (id >= nthreads)
+ return;
+
+ curandState rng = rng_states[id];
+ float3 poisiton = positions[id];
+ float3 direction = directions[id];
+ direction /= norm(direction);
+ float3 polarization = polarizations[id];
+ polarization /= norm(polarization);
+ float wavelength = wavelengths[id];
+ float last_hit_triangle = last_hit_triangles[id];
+
+ int steps = 0;
+ while (steps < max_steps)
+ {
+ steps++;
+
+ float distance;
+
+ last_hit_triangle = intersect_mesh(position, direction, start_node, first_node, distance, last_hit_triangle);
+
+ if (last_hit_triangle == -1)
+ break;
+
+ uint4 triangle_data = triangles[last_hit_triangle];
+
+ float3 v0 = vertices[triangle_data.x];
+ float3 v1 = vertices[triangle_data.y];
+ float3 v2 = vertices[triangle_data.z];
+
+ int material_in_index = convert(0xFF & (triangle_data.w >> 24));
+ int material_out_index = convert(0xFF & (triangle_data.w >> 16));
+ int surface_index = convert(0xFF & triangle_data.w >> 8);
+
+
+#endif
+
/* Trace the rays starting at `positions` traveling in the direction `directions`
to their intersection with the global mesh. If the ray intersects the mesh
set the pixel associated with the ray to a 32 bit color whose brightness is
@@ -370,22 +770,29 @@ __global__ void propagate(int nthreads, curandState *rng_states, float3 *positio
{
Surface surface = surfaces[surface_index];
- float absorption = interp_property(wavelength, surface.absorption);
- float reflection_diffuse = interp_property(wavelength, surface.reflection_diffuse);
- float reflection_specular = interp_property(wavelength, surface.reflection_specular);
+ float detect = interp_property(wavelength, surface.detect);
+ float absorb = interp_property(wavelength, surface.absorb);
+ float reflect_diffuse = interp_property(wavelength, surface.reflect_diffuse);
+ float reflect_specular = interp_property(wavelength, surface.reflect_specular);
// since the surface properties are interpolated linearly,
// we are guaranteed that they still sum to one.
float uniform_sample = curand_uniform(&rng);
- if (uniform_sample < absorption)
+ if (uniform_sample < absorb)
{
// absorb
state = SURFACE_ABSORB;
break;
}
- else if (uniform_sample < absorption + reflection_diffuse)
+ else if (uniform_sample < absorb + detect)
+ {
+ // detect
+ state = SURFACE_DETECT;
+ break;
+ }
+ else if (uniform_sample < absorb + detect + reflect_diffuse)
{
// diffusely reflect
direction = uniform_sphere(&rng);
@@ -465,4 +872,9 @@ __global__ void propagate(int nthreads, curandState *rng_states, float3 *positio
} // propagate
+#if 0
+
+
+#endif
+
} // extern "c"
diff --git a/src/materials.h b/src/materials.h
index 6115396..2355d24 100644
--- a/src/materials.h
+++ b/src/materials.h
@@ -15,9 +15,10 @@ struct Material
struct Surface
{
- float *absorption;
- float *reflection_diffuse;
- float *reflection_specular;
+ float *detect;
+ float *absorb;
+ float *reflect_diffuse;
+ float *reflect_specular;
};
__device__ Material materials[20];
@@ -54,11 +55,12 @@ __global__ void set_material(int material_index, float *refractive_index, float
materials[material_index].scattering_length = scattering_length;
}
-__global__ void set_surface(int surface_index, float *absorption, float *reflection_diffuse, float *reflection_specular)
+__global__ void set_surface(int surface_index, float *detect, float *absorb, float *reflect_diffuse, float *reflect_specular)
{
- surfaces[surface_index].absorption = absorption;
- surfaces[surface_index].reflection_diffuse = reflection_diffuse;
- surfaces[surface_index].reflection_specular = reflection_specular;
+ surfaces[surface_index].detect = detect;
+ surfaces[surface_index].absorb = absorb;
+ surfaces[surface_index].reflect_diffuse = reflect_diffuse;
+ surfaces[surface_index].reflect_specular = reflect_specular;
}
} // extern "c"
diff --git a/threadtest.py b/threadtest.py
index 92b0c81..2b72987 100644
--- a/threadtest.py
+++ b/threadtest.py
@@ -61,7 +61,7 @@ if __name__ == '__main__':
parser = optparse.OptionParser('%prog')
parser.add_option('-b', type='int', dest='nbits', default=8)
- parser.add_option('-j', type='string', dest='devices', default=1)
+ parser.add_option('-j', type='string', dest='devices', default='0')
parser.add_option('-n', type='int', dest='nblocks', default=64)
options, args = parser.parse_args()
diff --git a/view.py b/view.py
index 2abb76d..0a92dbd 100755
--- a/view.py
+++ b/view.py
@@ -2,6 +2,7 @@
import os
import sys
import numpy as np
+import pickle
import pygame
from pygame.locals import *
@@ -10,6 +11,7 @@ import src
from camera import *
from geometry import *
from transform import *
+#from bvh_view import *
import time
@@ -17,16 +19,7 @@ from pycuda import autoinit
from pycuda.compiler import SourceModule
from pycuda import gpuarray
-def view(viewable, name='', bits=8):
- """
- Render `viewable` in a pygame window.
-
- Movement:
- - zoom: scroll the mouse wheel
- - rotate: click and drag the mouse
- - move: shift+click and drag the mouse
- """
-
+def to_geometry(viewable, bits):
if isinstance(viewable, Geometry):
geometry = viewable
geometry.build(bits)
@@ -39,7 +32,81 @@ def view(viewable, name='', bits=8):
geometry.add_solid(Solid(viewable))
geometry.build(bits)
else:
- sys.exit("can't display %s" % args[0])
+ raise Exception("can't convert %s to a geometry" % viewable)
+
+ return geometry
+
+def box(lower_bound, upper_bound):
+ dx, dy, dz = upper_bound - lower_bound
+
+ vertices = np.array([lower_bound,
+ lower_bound + (dx,0,0),
+ lower_bound + (dx,dy,0),
+ lower_bound + (0,dy,0),
+ lower_bound + (0,0,dz),
+ lower_bound + (dx,0,dz),
+ lower_bound + (dx,dy,dz),
+ lower_bound + (0,dy,dz)])
+
+ triangles = np.empty((12,3), dtype=np.int32)
+ # bottom
+ triangles[0] = (1,3,2)
+ triangles[1] = (1,4,3)
+ # top
+ triangles[2] = (5,7,6)
+ triangles[3] = (5,8,7)
+ # left
+ triangles[4] = (5,1,2)
+ triangles[5] = (5,2,6)
+ # right
+ triangles[6] = (3,4,8)
+ triangles[7] = (3,8,7)
+ # front
+ triangles[8] = (2,3,7)
+ triangles[9] = (2,7,6)
+ # back
+ triangles[10] = (1,5,8)
+ triangles[11] = (1,8,4)
+
+ triangles -= 1
+
+ return Mesh(vertices, triangles)
+
+def bvh_mesh(geometry, layer):
+ lower_bounds = geometry.lower_bounds[geometry.layers == layer]
+ upper_bounds = geometry.upper_bounds[geometry.layers == layer]
+
+ if len(lower_bounds) == 0 or len(upper_bounds) == 0:
+ raise Exception('no nodes at layer %i' % layer)
+
+ mesh = box(lower_bounds[0], upper_bounds[0])
+
+ for lower_bound, upper_bound in zip(lower_bounds, upper_bounds)[1:]:
+ mesh += box(lower_bound, upper_bound)
+
+ return mesh
+
+def view(viewable, name='', bits=8, load_bvh=False, camera_view=None):
+ """
+ Render `viewable` in a pygame window.
+
+ Movement:
+ - zoom: scroll the mouse wheel
+ - rotate: click and drag the mouse
+ - move: shift+click and drag the mouse
+ """
+
+ geometry = to_geometry(viewable, bits)
+
+ if load_bvh:
+ print 'loading bvh...'
+ bvhg = []
+
+ bvhg.append(geometry)
+
+ for layer in sorted(np.unique(geometry.layers)):
+ print 'building layer %i' % layer
+ bvhg.append(to_geometry(bvh_mesh(geometry, layer), bits))
lower_bound = np.array([np.min(geometry.mesh[:][:,:,0]),
np.min(geometry.mesh[:][:,:,1]),
@@ -65,27 +132,36 @@ def view(viewable, name='', bits=8):
screen = pygame.display.set_mode(size)
pygame.display.set_caption(name)
- camera = Camera(size)
+ if camera_view is None:
+ camera = Camera(size)
- diagonal = np.linalg.norm(upper_bound-lower_bound)
+ diagonal = np.linalg.norm(upper_bound-lower_bound)
- 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)
+ 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 = camera.get_rays()
- origins_float3 = np.empty(origins.shape[0], dtype=gpuarray.vec.float3)
- origins_float3['x'] = origins[:,0]
- origins_float3['y'] = origins[:,1]
- origins_float3['z'] = origins[:,2]
+ origins_float3 = np.empty(origins.shape[0], dtype=gpuarray.vec.float3)
+ origins_float3['x'] = origins[:,0]
+ origins_float3['y'] = origins[:,1]
+ origins_float3['z'] = origins[:,2]
- directions_float3 = np.empty(directions.shape[0], dtype=gpuarray.vec.float3)
- directions_float3['x'] = directions[:,0]
- directions_float3['y'] = directions[:,1]
- directions_float3['z'] = directions[:,2]
+ directions_float3 = np.empty(directions.shape[0], dtype=gpuarray.vec.float3)
+ directions_float3['x'] = directions[:,0]
+ directions_float3['y'] = directions[:,1]
+ directions_float3['z'] = directions[:,2]
+ else:
+ f = open(camera_view, 'rb')
+ origins_float3 = pickle.load(f)
+ directions_float3 = pickle.load(f)
+ point = pickle.load(f)
+ axis1 = pickle.load(f)
+ axis2 = pickle.load(f)
+ f.close()
origins_gpu = cuda.to_device(origins_float3)
directions_gpu = cuda.to_device(directions_float3)
@@ -104,7 +180,7 @@ def view(viewable, name='', bits=8):
cuda.Context.synchronize()
elapsed = time.time() - t0
- print 'elapsed %f sec' % elapsed
+ #print 'elapsed %f sec' % elapsed
cuda.memcpy_dtoh(pixels, pixels_gpu)
pygame.surfarray.blit_array(screen, pixels.reshape(size))
@@ -116,6 +192,8 @@ def view(viewable, name='', bits=8):
clicked = False
shift = False
+ current_layer = 0
+
while not done:
for event in pygame.event.get():
if event.type == MOUSEBUTTONDOWN:
@@ -187,6 +265,57 @@ def view(viewable, name='', bits=8):
done = True
break
+ if event.key == K_PAGEUP:
+ try:
+ if current_layer+1 >= len(bvhg):
+ raise IndexError
+
+ geometry = bvhg[current_layer+1]
+ current_layer += 1
+
+ geometry.load(module, color=True)
+ render()
+ except IndexError:
+ print 'no further layers to view'
+
+ if event.key == K_PAGEDOWN:
+ try:
+ if current_layer-1 < 0:
+ raise IndexError
+
+ geometry = bvhg[current_layer-1]
+ current_layer -= 1
+
+ geometry.load(module, color=True)
+ render()
+ except IndexError:
+ print 'no further layers to view'
+
+ if event.key == K_F11:
+ if name == '':
+ root, ext = 'camera_view', 'pkl'
+ else:
+ root, ext = name + '_camera_view', 'pkl'
+
+ filename = '.'.join([root, ext])
+
+ i = 1
+ while os.path.exists(filename):
+ filename = '.'.join([root + str(i), ext])
+ i += 1
+
+ f = open(filename, 'wb')
+ cuda.memcpy_dtoh(origins_float3, origins_gpu)
+ cuda.memcpy_dtoh(directions_float3, directions_gpu)
+ pickle.dump(origins_float3, f)
+ pickle.dump(directions_float3, f)
+ pickle.dump(point, f)
+ pickle.dump(axis1, f)
+ pickle.dump(axis2, f)
+ f.close()
+
+ print 'camera view saved to %s' % filename
+
if event.key == K_F12:
if name == '':
root, ext = 'screenshot', 'png'
@@ -220,6 +349,10 @@ if __name__ == '__main__':
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('-v', '--camera_view', dest='camera_view',
+ help='load a camera view', default=None)
options, args = parser.parse_args()
if len(args) < 1:
@@ -229,16 +362,13 @@ if __name__ == '__main__':
root, ext = os.path.splitext(tail)
if ext.lower() == '.stl':
- geometry = Geometry()
- geometry.add_solid(Solid(mesh_from_stl(args[0])))
- geometry.build(options.bits)
- view(geometry, tail)
+ view(mesh_from_stl(args[0]), root, options.bits, options.load_bvh, options.camera_view)
else:
import inspect
members = dict(inspect.getmembers(detectors) + inspect.getmembers(solids))
if args[0] in members:
- view(members[args[0]], args[0], options.bits)
+ view(members[args[0]], args[0], options.bits, options.load_bvh, options.camera_view)
else:
sys.exit("couldn't find object %s" % args[0])
generated by cgit (git 2.34.1) at 2025-05-08 17:24:25 +0000