added a hybrid monte carlo ray tracing rendering algorithm

1 files changed, 12 insertions, 15 deletions

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))