fix nan bug by clamping dot() of surface normal and -photon direction to [0.0,1.0]

1 files changed, 7 insertions, 2 deletions

diff --git a/src/photon.h b/src/photon.h
index 4203a50..6cf3367 100644
--- a/src/photon.h
+++ b/src/photon.h
@@ -50,6 +50,11 @@ enum
 
 enum {BREAK, CONTINUE, PASS}; // return value from propagate_to_boundary
 
+__device__ float get_theta(const float3 &a, const float3 &b)
+{
+	return acosf(fmaxf(0.0f,fminf(1.0f,dot(a,b))));
+}
+
 __device__ void fill_state(State &s, Photon &p)
 {
 	p.last_hit_triangle = intersect_mesh(p.position, p.direction, s.distance_to_boundary, p.last_hit_triangle);
@@ -169,7 +174,7 @@ __device__ int propagate_to_boundary(Photon &p, State &s, curandState &rng)
 
 __device__ void propagate_at_boundary(Photon &p, State &s, curandState &rng)
 {
-	float incident_angle = acosf(dot(s.surface_normal, -p.direction));
+	float incident_angle = get_theta(s.surface_normal,-p.direction);
 	float refracted_angle = asinf(sinf(incident_angle)*s.refractive_index1/s.refractive_index2);
 
 	float3 incident_plane_normal = cross(p.direction, s.surface_normal);
@@ -260,7 +265,7 @@ __device__ int propagate_at_surface(Photon &p, State &s, curandState &rng)
 	else
 	{
 		// specularly reflect
-		float incident_angle = acosf(dot(s.surface_normal, -p.direction));
+		float incident_angle = get_theta(s.surface_normal,-p.direction);
 		float3 incident_plane_normal = cross(p.direction, s.surface_normal);
 		incident_plane_normal /= norm(incident_plane_normal);