simplify surface models

remove the ``SURFACE_SPECULAR`` and ``SURFACE_DIFFUSE`` models, since their functionality is available using the more-general ``SURFACE_DEFAULT``. also allow the user to specify the reflection type (specular/diffuse) for the complex and wls models. change wls so the normalization of properties is more consistent with the default.

4 files changed, 33 insertions, 32 deletions

diff --git a/chroma/cuda/geometry_types.h b/chroma/cuda/geometry_types.h
index f82f372..46226b2 100644
--- a/chroma/cuda/geometry_types.h
+++ b/chroma/cuda/geometry_types.h
@@ -14,8 +14,6 @@ struct Material
 // surface models
 enum {
     SURFACE_DEFAULT,  // specular + diffuse + absorption + detection
-    SURFACE_SPECULAR, // perfect specular reflector
-    SURFACE_DIFFUSE,  // perfect diffuse reflector
     SURFACE_COMPLEX,  // use complex index of refraction
     SURFACE_WLS       // wavelength-shifting reemission
 };
@@ -25,7 +23,6 @@ struct Surface
     float *detect;
     float *absorb;
     float *reemit;
-    float *reflect;
     float *reflect_diffuse;
     float *reflect_specular;
     float *eta;
@@ -71,3 +68,4 @@ struct Geometry
 };
 
 #endif
+
diff --git a/chroma/cuda/photon.h b/chroma/cuda/photon.h
index 09430e9..4ec1cea 100644
--- a/chroma/cuda/photon.h
+++ b/chroma/cuda/photon.h
@@ -332,7 +332,7 @@ propagate_at_specular_reflector(Photon &p, State &s)
     p.history |= REFLECT_SPECULAR;
 
     return CONTINUE;
-} // reflect_specular
+} // propagate_at_specular_reflector
 
 __device__ int
 propagate_at_diffuse_reflector(Photon &p, State &s, curandState &rng)
@@ -348,17 +348,17 @@ propagate_at_diffuse_reflector(Photon &p, State &s, curandState &rng)
     p.history |= REFLECT_DIFFUSE;
 
     return CONTINUE;
-} // reflect_diffuse
+} // propagate_at_diffuse_reflector
 
 __device__ int
 propagate_complex(Photon &p, State &s, curandState &rng, Surface* surface, bool use_weights=false)
 {
     float detect = interp_property(surface, p.wavelength, surface->detect);
+    float reflect_specular = interp_property(surface, p.wavelength, surface->reflect_specular);
+    float reflect_diffuse = interp_property(surface, p.wavelength, surface->reflect_diffuse);
     float n2_eta = interp_property(surface, p.wavelength, surface->eta);
     float n2_k = interp_property(surface, p.wavelength, surface->k);
 
-    float uniform_sample = curand_uniform(&rng);
-
     // thin film optical model, adapted from RAT PMT optical model by P. Jones
     cuFloatComplex n1 = make_cuFloatComplex(s.refractive_index1, 0.0f);
     cuFloatComplex n2 = make_cuFloatComplex(n2_eta, n2_k);
@@ -475,28 +475,33 @@ propagate_complex(Photon &p, State &s, curandState &rng, Surface* surface, bool
     }
 
     if (use_weights && detect > 0.0f) {
-        p.history |= SURFACE_ABSORB;
         p.history |= SURFACE_DETECT;
         p.weight *= detect;
         return BREAK;
     }
 
-    if (uniform_sample < absorb) {
-        // absorb
-        p.history |= SURFACE_ABSORB;
+    float uniform_sample = curand_uniform(&rng);
 
+    if (uniform_sample < absorb) {
         // detection probability is conditional on absorption here
         float uniform_sample_detect = curand_uniform(&rng);
         if (uniform_sample_detect < detect)
             p.history |= SURFACE_DETECT;
+        else
+            p.history |= SURFACE_ABSORB;
 
         return BREAK;
     }
-    else if (uniform_sample < absorb + reflect) {
-        return propagate_at_diffuse_reflector(p, s, rng);
+    else if (uniform_sample < absorb + reflect || !surface->transmissive) {
+        // reflect, specularly (default) or diffusely
+        float uniform_sample_reflect = curand_uniform(&rng);
+        if (uniform_sample_reflect < reflect_diffuse)
+            return propagate_at_diffuse_reflector(p, s, rng);
+        else
+            return propagate_at_specular_reflector(p, s);
     }
     else {
-        // transmit
+        // refract and transmit
         p.direction = rotate(s.surface_normal, PI-refracted_angle, incident_plane_normal);
         p.polarization = cross(incident_plane_normal, p.direction);
         p.polarization /= norm(p.polarization);
@@ -509,7 +514,8 @@ __device__ int
 propagate_at_wls(Photon &p, State &s, curandState &rng, Surface *surface, bool use_weights=false)
 {
     float absorb = interp_property(surface, p.wavelength, surface->absorb);
-    float reflect = interp_property(surface, p.wavelength, surface->reflect);
+    float reflect_specular = interp_property(surface, p.wavelength, surface->reflect_specular);
+    float reflect_diffuse = interp_property(surface, p.wavelength, surface->reflect_diffuse);
     float reemit = interp_property(surface, p.wavelength, surface->reemit);
 
     float uniform_sample = curand_uniform(&rng);
@@ -519,12 +525,11 @@ propagate_at_wls(Photon &p, State &s, curandState &rng, Surface *surface, bool u
         float survive = 1.0f - absorb;
         absorb = 0.0f;
         p.weight *= survive;
-        reflect /= survive;
+        reflect_diffuse /= survive;
+        reflect_specular /= survive;
     }
 
     if (uniform_sample < absorb) {
-        p.history |= SURFACE_ABSORB;
-
         float uniform_sample_reemit = curand_uniform(&rng);
         if (uniform_sample_reemit < reemit) {
             p.history |= SURFACE_REEMIT;
@@ -532,10 +537,16 @@ propagate_at_wls(Photon &p, State &s, curandState &rng, Surface *surface, bool u
             return propagate_at_diffuse_reflector(p, s, rng); // reemit isotropically (eh?)
         }
 
+        p.history |= SURFACE_ABSORB;
         return BREAK;
     }
-    else if (uniform_sample < absorb + reflect) {
-        return propagate_at_diffuse_reflector(p, s, rng);
+    else if (uniform_sample < absorb + reflect_specular + reflect_diffuse) {
+        // choose how to reflect, defaulting to diffuse
+        float uniform_sample_reflect = curand_uniform(&rng) * (reflect_specular + reflect_diffuse);
+        if (uniform_sample_reflect < reflect_specular)
+            return propagate_at_specular_reflector(p, s);
+        else
+            return propagate_at_diffuse_reflector(p, s, rng);
     }
     else {
         p.history |= SURFACE_TRANSMIT;
@@ -549,18 +560,14 @@ propagate_at_surface(Photon &p, State &s, curandState &rng, Geometry *geometry,
 {
     Surface *surface = geometry->surfaces[s.surface_index];
 
-    if (surface->model == SURFACE_SPECULAR)
-        return propagate_at_specular_reflector(p, s);
-    else if (surface->model == SURFACE_DIFFUSE)
-        return propagate_at_diffuse_reflector(p, s, rng);
-    else if (surface->model == SURFACE_COMPLEX)
+    if (surface->model == SURFACE_COMPLEX)
         return propagate_complex(p, s, rng, surface, use_weights);
     else if (surface->model == SURFACE_WLS)
         return propagate_at_wls(p, s, rng, surface, use_weights);
     else {
-        // if no surface model is specified, do a combination of specular and
+        // use default surface model: do a combination of specular and
         // diffuse reflection, detection, and absorption based on relative
-        // probabilties (i.e. the "old" behavior)
+        // probabilties
 
         // since the surface properties are interpolated linearly, we are
         // guaranteed that they still sum to 1.0.
diff --git a/chroma/geometry.py b/chroma/geometry.py
index 416a6e7..19294be 100644
--- a/chroma/geometry.py
+++ b/chroma/geometry.py
@@ -175,7 +175,6 @@ class Surface(object):
         self.set('detect', 0)
         self.set('absorb', 0)
         self.set('reemit', 0)
-        self.set('reflect', 0)
         self.set('reflect_diffuse', 0)
         self.set('reflect_specular', 0)
         self.set('eta', 0)
diff --git a/chroma/gpu/geometry.py b/chroma/gpu/geometry.py
index f74d32d..ff71b56 100644
--- a/chroma/gpu/geometry.py
+++ b/chroma/gpu/geometry.py
@@ -83,8 +83,6 @@ class GPUGeometry(object):
             absorb_gpu = ga.to_gpu(absorb)
             reemit = interp_material_property(wavelengths, surface.reemit)
             reemit_gpu = ga.to_gpu(reemit)
-            reflect = interp_material_property(wavelengths, surface.reflect)
-            reflect_gpu = ga.to_gpu(reflect)
             reflect_diffuse = interp_material_property(wavelengths, surface.reflect_diffuse)
             reflect_diffuse_gpu = ga.to_gpu(reflect_diffuse)
             reflect_specular = interp_material_property(wavelengths, surface.reflect_specular)
@@ -99,7 +97,6 @@ class GPUGeometry(object):
             self.surface_data.append(detect_gpu)
             self.surface_data.append(absorb_gpu)
             self.surface_data.append(reemit_gpu)
-            self.surface_data.append(reflect_gpu)
             self.surface_data.append(reflect_diffuse_gpu)
             self.surface_data.append(reflect_specular_gpu)
             self.surface_data.append(eta_gpu)
@@ -108,7 +105,7 @@ class GPUGeometry(object):
 
             surface_gpu = \
                 make_gpu_struct(surface_struct_size,
-                                [detect_gpu, absorb_gpu, reemit_gpu, reflect_gpu,
+                                [detect_gpu, absorb_gpu, reemit_gpu,
                                  reflect_diffuse_gpu,reflect_specular_gpu,
                                  eta_gpu, k_gpu, reemission_cdf_gpu,
                                  np.uint32(surface.model),