speedup mesh intersection by skipping directly to the child node of nodes with only one child.

3 files changed, 27 insertions, 25 deletions

diff --git a/geometry.py b/geometry.py
index 32d6462..2405e07 100644
--- a/geometry.py
+++ b/geometry.py
@@ -283,7 +283,7 @@ class Geometry(object):
         self.first_node = unique_zvalues.size
 
         begin_last_layer = 0
-        
+
         while True:
             bit_shifted_zvalues = zvalues >> 1
             unique_bit_shifted_zvalues = np.unique(bit_shifted_zvalues)
diff --git a/src/kernel.cu b/src/kernel.cu
index 55fde70..f257c6c 100644
--- a/src/kernel.cu
+++ b/src/kernel.cu
@@ -81,34 +81,37 @@ __device__ int intersect_mesh(const float3 &origin, const float3& direction, con
 
 	do
 	{
-		int index = tex1Dfetch(node_map, *node);
-		int length = tex1Dfetch(node_length, *node);
+		int first_child = tex1Dfetch(node_map, *node);
+		int child_count = tex1Dfetch(node_length, *node);
 
+		while (*node >= first_node && child_count == 1)
+		{
+			*node = first_child;
+			first_child = tex1Dfetch(node_map, *node);
+			child_count = tex1Dfetch(node_length, *node);
+		}
+		
 		if (*node >= first_node)
 		{
-			for (i=0; i < length; i++)
+			for (i=0; i < child_count; i++)
 			{
-				if (intersect_node(origin, direction, index+i))
+				if (intersect_node(origin, direction, first_child+i))
 				{
-					//*node++ = index+i;
-					*node = index+i;
+					*node = first_child+i;
 					node++;
 				}
 			}
 
-			if (node == head)
-				break;
-
 			node--;
 		}
 		else // node is a leaf
 		{
-			for (i=0; i < length; i++)
+			for (i=0; i < child_count; i++)
 			{
-				if (last_hit_triangle == index+i)
+				if (last_hit_triangle == first_child+i)
 					continue;
 
-				uint4 triangle_data = triangles[index+i];
+				uint4 triangle_data = triangles[first_child+i];
 
 				float3 v0 = vertices[triangle_data.x];
 				float3 v1 = vertices[triangle_data.y];
@@ -118,14 +121,14 @@ __device__ int intersect_mesh(const float3 &origin, const float3& direction, con
 				{
 					if (triangle_index == -1)
 					{
-						triangle_index = index + i;
+						triangle_index = first_child + i;
 						min_distance = distance;
 						continue;
 					}
 
 					if (distance < min_distance)
 					{
-						triangle_index = index + i;
+						triangle_index = first_child + i;
 						min_distance = distance;
 					}
 				}
@@ -245,7 +248,9 @@ __global__ void propagate(int nthreads, curandState *rng_states, float3 *positio
 	curandState rng = rng_states[id];
 	float3 position = positions[id];
 	float3 direction = directions[id];
+	direction /= norm(direction);
 	float3 polarization = polarizations[id];
+	polarization /= norm(polarization);
 	float wavelength = wavelengths[id];
 	float time = times[id];
 	int last_hit_triangle = last_hit_triangles[id];
@@ -255,9 +260,6 @@ __global__ void propagate(int nthreads, curandState *rng_states, float3 *positio
 	{
 		steps++;
 
-		direction /= norm(direction);
-		polarization /= norm(polarization);
-
 		float distance;
 
 	        last_hit_triangle = intersect_mesh(position, direction, start_node, first_node, distance, last_hit_triangle);
@@ -372,10 +374,8 @@ __global__ void propagate(int nthreads, curandState *rng_states, float3 *positio
 			float reflection_diffuse = interp_property(wavelength, surface.reflection_diffuse);
 			float reflection_specular = interp_property(wavelength, surface.reflection_specular);
 
-			float sum = absorption + reflection_diffuse + reflection_specular;
-			absorption /= sum;
-			reflection_diffuse /= sum;
-			reflection_specular /= sum;
+			// since the surface properties are interpolated linearly,
+			// we are guaranteed that they still sum to one.
 
 			float uniform_sample = curand_uniform(&rng);
 
@@ -402,9 +402,10 @@ __global__ void propagate(int nthreads, curandState *rng_states, float3 *positio
 			else
 			{
 				// specularly reflect
-
 				direction = rotate(surface_normal, incident_angle, p);
 
+				// polarization ?
+
 				continue;
 			}
 
@@ -446,6 +447,7 @@ __global__ void propagate(int nthreads, curandState *rng_states, float3 *positio
 			}
 
 			polarization = cross(p, direction);
+			polarization /= norm(polarization);
 
 			continue;
 		} // photon polarization parallel to surface
diff --git a/threadtest.py b/threadtest.py
index 6cc9c07..5dffb98 100644
--- a/threadtest.py
+++ b/threadtest.py
@@ -35,7 +35,7 @@ def likelihood(detector, event_times, position=(0,0,0), nphotons=5000, neval=100
     log_likelihood_variance = 0.0
     for i, time in event_times:
         h = Histogram(500, (-0.5e-9, 99.5e-9))
-        h.fill(t[:,i][t[:,i] < 1e8])
+        h.fill(t[:,i])
 
         if h.nentries > 0:
             h.normalize()
@@ -67,7 +67,7 @@ if __name__ == '__main__':
 
     from detectors import lbne_cut, minilbne_cut, minilbne, microlbne
 
-    detector = minilbne_cut
+    detector = minilbne
 
     detector.build(bits=options.nbits)