fix(CSGOptiX): handle coincident sibling faces

CSGOptiX/CMakeLists.txt

@@ -22,6 +22,7 @@ list(APPEND SOURCES
 list(APPEND HEADERS
     CSGOptiX.h
     CSGOptiXService.h
+    CSGOptiXHelpers.h
     CSGOPTIX_API_EXPORT.hh
     CSGOPTIX_LOG.hh
 )

CSGOptiX/CSGOptiX.cc

@@ -592,25 +592,14 @@ void CSGOptiX::initSimulate()
     params->PropagateEpsilon0Mask = SEventConfig::PropagateEpsilon0Mask();  // eg from CK|SI|TO|SC|RE
 
     params->PropagateRefine = SEventConfig::PropagateRefine();
-    params->PropagateRefineDistance = SEventConfig::PropagateRefineDistance();  // approx distance beyond which to refine intersect with 2nd trace
-
-    params->tmin = SEventConfig::PropagateEpsilon() ;  // eg 0.1 0.05 to avoid self-intersection off boundaries
-    params->tmax = 1000000.f ;
-    params->max_time = SEventConfig::MaxTime() ;
+    params->PropagateRefineDistance = SEventConfig::PropagateRefineDistance(); // approx distance beyond which to refine intersect with 2nd trace
 
+    params->tmin = SEventConfig::PropagateEpsilon(); // eg 0.1 0.05 to avoid self-intersection off boundaries
 
+    params->tmax = 1000000.f;
+    params->max_time = SEventConfig::MaxTime();
 }
 
-
-
-
-
-
-
-
-
-
-
 /**
 CSGOptiX::initRender
 --------------------------

CSGOptiX/CSGOptiX7.cu

@@ -77,6 +77,7 @@ __intersection__is
 
 #include "Binding.h"
 #include "Params.h"
+#include "CSGOptiXHelpers.h"
 
 #ifdef WITH_PRD
 #include "scuda_pointer.h"
@@ -439,6 +440,12 @@ static __forceinline__ __device__ void simulate( const uint3& launch_idx, const
 #endif
 
     sim->generate_photon(ctx.p, rng, gs, photon_idx, genstep_idx );
+    // Only Cerenkov and scintillation gensteps store a source-material line in
+    // q0.u.z. TORCH and FRAME use that slot for a genstep id, so leave the
+    // carried line invalid for all other genstep types.
+    const unsigned gentype = gs.q0.u.x;
+    const bool carries_matline = CSGOptiX7_GenstepCarriesMaterialLine(gentype);
+    ctx.current_matline = carries_matline ? gs.q0.u.z : 0xFFFFFFFFu;
 
     FlowAction command = FlowAction::Start;
     int bounce = 0 ;
@@ -902,8 +909,16 @@ extern "C" __global__ void __intersection__is()
         const unsigned boundary = node->boundary() ;  // all CSGNode in the tree for one CSGPrim tree have same boundary
         const unsigned globalPrimIdx_boundary = (( globalPrimIdx & 0xffffu ) << 16 ) | ( boundary & 0xffffu ) ;
 
+        // Coincident sibling faces can produce equal entering and exiting hits.
+        // Move only the reported exiting distance to the next representable float
+        // during photon simulation so the entering sibling wins an exact tie.
+        // This one-ULP ordering change cannot overstep a representable physical
+        // gap. Keep the unbiased isect.w in the PRD for position advancement.
+        const float cosI = dot(ray_direction, make_float3(isect.x, isect.y, isect.z));
+        const float t_report = CSGOptiX7_ReportedIntersectionDistance(isect.w, params.raygenmode, cosI);
+
 #ifdef WITH_PRD
-        if(optixReportIntersection( isect.w, hitKind))
+        if (optixReportIntersection(t_report, hitKind))
         {
             quad2* prd = SOPTIX_getPRD<quad2>(); // access prd addr from RG program
             prd->q0.f = isect ;  // .w:distance and .xyz:normal which starts as the local frame one
@@ -919,7 +934,7 @@ extern "C" __global__ void __intersection__is()
         a3 = __float_as_uint( isect.w ) ;
         a4 = globalPrimIdx_boundary ;
         a5 = __float_as_uint( lposcost );
-        optixReportIntersection( isect.w, hitKind, a0, a1, a2, a3, a4, a5 );
+        optixReportIntersection(t_report, hitKind, a0, a1, a2, a3, a4, a5);
 
         // IS:optixReportIntersection writes the attributes that can be read in CH and AH programs
         // max 8 attribute registers, see PIP::PIP, communicate to __closesthit__ch

CSGOptiX/CSGOptiXHelpers.h

@@ -0,0 +1,30 @@
+#pragma once
+
+#include <math.h>
+
+#include "OpticksGenstep.h"
+#include "SRG.h"
+
+#if defined(__CUDACC__) || defined(__CUDABE__)
+#define CSGOPTIX7_HELPER_METHOD __device__ __forceinline__
+#else
+#define CSGOPTIX7_HELPER_METHOD inline
+#endif
+
+CSGOPTIX7_HELPER_METHOD bool CSGOptiX7_GenstepCarriesMaterialLine(unsigned gentype)
+{
+    return gentype == OpticksGenstep_CERENKOV ||
+           gentype == OpticksGenstep_SCINTILLATION ||
+           gentype == OpticksGenstep_G4Cerenkov_modified;
+}
+
+CSGOPTIX7_HELPER_METHOD float CSGOptiX7_ReportedIntersectionDistance(
+    float    distance,
+    unsigned raygenmode,
+    float    direction_dot_normal)
+{
+    const bool defer_exit = raygenmode == SRG_SIMULATE && direction_dot_normal > 0.f;
+    return defer_exit ? nextafterf(distance, INFINITY) : distance;
+}
+
+#undef CSGOPTIX7_HELPER_METHOD

CSGOptiX/Params.cc

@@ -123,12 +123,10 @@ std::string Params::detail() const
         << std::setw(20) << " evt " << std::setw(10) << evt  << std::endl
         ;
     std::string s = ss.str();
-    return s ;
+    return s;
 }
 
-
-Params::Params(int raygenmode_, unsigned width, unsigned height, unsigned depth)
-    :
+Params::Params(int raygenmode_, unsigned width, unsigned height, unsigned depth) :
     raygenmode(SRG_RENDER),
     node(nullptr),
     plan(nullptr),

CSGOptiX/Params.h

@@ -61,12 +61,12 @@ struct Params
     float3     V ;
     float3     W ;
     float3     WNORM ;
-    float      tmin ;
+    float      tmin;
 
-    float      tmin0 ;
-    unsigned   PropagateEpsilon0Mask ;   // default from SEventConfig TO,CK,SI,SC,RE
-    float      PropagateRefineDistance ;
-    bool       PropagateRefine ;
+    float    tmin0;
+    unsigned PropagateEpsilon0Mask; // default from SEventConfig TO,CK,SI,SC,RE
+    float    PropagateRefineDistance;
+    bool     PropagateRefine;
 
     float      tmax ;
     float4     ZPROJ ;

CSGOptiX/tests/CMakeLists.txt

@@ -10,6 +10,7 @@ or doing geometry specific things.
 
 set(SINGLE_SOURCES
 
+    CSGOptiXHelpersTest.cu
     CSGOptiXDescTest.cc
     CSGOptiXSimtraceTest.cc
     CSGOptiXSimulateTest.cc
@@ -49,4 +50,6 @@ foreach(SRC ${TEST_SOURCES})
     )
 endforeach()
 
+add_test(NAME ${name}.CSGOptiXHelpersTest COMMAND $<TARGET_FILE:CSGOptiXHelpersTest>)
+
 set_tests_properties(${name}.CSGOptiXRenderTest PROPERTIES DISABLED TRUE)

CSGOptiX/tests/CSGOptiXHelpersTest.cu

@@ -0,0 +1,92 @@
+#include <cstdlib>
+#include <cstdio>
+
+#include <cuda_runtime.h>
+
+#include "CSGOptiXHelpers.h"
+
+namespace
+{
+    void require(bool condition, const char* expression, int line)
+    {
+        if (condition) return;
+        std::fprintf(stderr, "%s:%d requirement failed: %s\n", __FILE__, line, expression);
+        std::exit(EXIT_FAILURE);
+    }
+}
+
+#define REQUIRE(expression) require((expression), #expression, __LINE__)
+
+struct CSGOptiXHelpersTestResult
+{
+    bool torch_carries_matline;
+    bool frame_carries_matline;
+    bool cerenkov_carries_matline;
+    bool scintillation_carries_matline;
+    bool modified_cerenkov_carries_matline;
+    float entering_distance;
+    float exiting_distance;
+    float non_simulation_exit_distance;
+};
+
+__global__ void CSGOptiXHelpersTest_kernel(CSGOptiXHelpersTestResult* result)
+{
+    result->torch_carries_matline =
+        CSGOptiX7_GenstepCarriesMaterialLine(OpticksGenstep_TORCH);
+    result->frame_carries_matline =
+        CSGOptiX7_GenstepCarriesMaterialLine(OpticksGenstep_FRAME);
+    result->cerenkov_carries_matline =
+        CSGOptiX7_GenstepCarriesMaterialLine(OpticksGenstep_CERENKOV);
+    result->scintillation_carries_matline =
+        CSGOptiX7_GenstepCarriesMaterialLine(OpticksGenstep_SCINTILLATION);
+    result->modified_cerenkov_carries_matline =
+        CSGOptiX7_GenstepCarriesMaterialLine(OpticksGenstep_G4Cerenkov_modified);
+
+    const float coincident_distance = 10.f;
+    result->entering_distance = CSGOptiX7_ReportedIntersectionDistance(
+        coincident_distance,
+        SRG_SIMULATE,
+        -1.f);
+    result->exiting_distance = CSGOptiX7_ReportedIntersectionDistance(
+        coincident_distance,
+        SRG_SIMULATE,
+        1.f);
+    result->non_simulation_exit_distance = CSGOptiX7_ReportedIntersectionDistance(
+        coincident_distance,
+        SRG_SIMTRACE,
+        1.f);
+}
+
+int main()
+{
+    CSGOptiXHelpersTestResult* device_result = nullptr;
+    REQUIRE(cudaMalloc(&device_result, sizeof(CSGOptiXHelpersTestResult)) == cudaSuccess);
+
+    CSGOptiXHelpersTest_kernel<<<1, 1>>>(device_result);
+    REQUIRE(cudaGetLastError() == cudaSuccess);
+    REQUIRE(cudaDeviceSynchronize() == cudaSuccess);
+
+    CSGOptiXHelpersTestResult result = {};
+    REQUIRE(cudaMemcpy(
+        &result,
+        device_result,
+        sizeof(CSGOptiXHelpersTestResult),
+        cudaMemcpyDeviceToHost) == cudaSuccess);
+    REQUIRE(cudaFree(device_result) == cudaSuccess);
+
+    REQUIRE(!result.torch_carries_matline);
+    REQUIRE(!result.frame_carries_matline);
+    REQUIRE(result.cerenkov_carries_matline);
+    REQUIRE(result.scintillation_carries_matline);
+    REQUIRE(result.modified_cerenkov_carries_matline);
+
+    REQUIRE(result.entering_distance == 10.f);
+    REQUIRE(result.exiting_distance == nextafterf(10.f, INFINITY));
+    REQUIRE(result.entering_distance < result.exiting_distance);
+    REQUIRE(result.non_simulation_exit_distance == 10.f);
+
+    std::puts("CSGOptiXHelpersTest: PASS");
+    return 0;
+}
+
+#undef REQUIRE

config/sibling_pair.json

@@ -0,0 +1,34 @@
+{
+  "torch": {
+    "gentype": "TORCH",
+    "trackid": 0,
+    "matline": 0,
+    "numphoton": 2000,
+
+    "pos": [-10.0, 0.0, -50.1],
+    "time": 0.0,
+
+    "mom": [0.3, 0.0, 1.0],
+    "weight": 0.0,
+
+    "pol": [0.0, 1.0, 0.0],
+    "wavelength": 420.0,
+
+    "zenith": [-1.0, 1.0],
+    "azimuth": [0.25, 0.25],
+
+    "radius": 0.05,
+    "distance": 0.0,
+    "mode": 255,
+    "type": "disc"
+  },
+
+  "event": {
+    "mode": "HitPhoton",
+    "maxslot": 10000,
+    "max_bounce": 8,
+    "propagate_epsilon": 0.05,
+    "propagate_epsilon0": 0.05,
+    "propagate_epsilon0_mask": "TO,CK,SI,SC,RE"
+  }
+}

qudarap/qbnd.h

@@ -35,16 +35,13 @@ struct qbnd
     unsigned            boundary_tex_MaterialLine_LS ;
     quad*               optical ;
 
-#if defined(__CUDACC__) || defined(__CUDABE__) || defined( MOCK_TEXTURE) || defined(MOCK_CUDA)
-    QBND_METHOD float4  boundary_lookup( unsigned ix, unsigned iy );
-    QBND_METHOD float4  boundary_lookup( float nm, unsigned line, unsigned k );
-    QBND_METHOD void    fill_state(sstate& s, unsigned boundary, float wavelength, float cosTheta, unsigned long long idx, unsigned long long base_pidx );
+#if defined(__CUDACC__) || defined(__CUDABE__) || defined(MOCK_TEXTURE) || defined(MOCK_CUDA)
+    QBND_METHOD float4 boundary_lookup(unsigned ix, unsigned iy);
+    QBND_METHOD float4 boundary_lookup(float nm, unsigned line, unsigned k);
+    QBND_METHOD void fill_state(sstate& s, unsigned boundary, float wavelength, float cosTheta, unsigned long long idx, unsigned long long base_pidx, unsigned carried_matline = 0xFFFFFFFFu);
 #endif
-
 };
 
-
-
 #if defined( MOCK_TEXTURE) || defined(MOCK_CUDA)
 #include "stexture.h"
 #endif
@@ -194,30 +191,51 @@ s.optical.x
 
 **/
 
-inline QBND_METHOD void qbnd::fill_state(sstate& s, unsigned boundary, float wavelength, float cosTheta, unsigned long long idx, unsigned long long base_pidx  )
+inline QBND_METHOD void qbnd::fill_state(sstate& s, unsigned boundary, float wavelength, float cosTheta, unsigned long long idx, unsigned long long base_pidx, unsigned carried_matline)
 {
-    const int line = boundary*_BOUNDARY_NUM_MATSUR ;      // now that are not signing boundary use 0-based
-
-    const int m1_line = cosTheta > 0.f ? line + IMAT : line + OMAT ;
-    const int m2_line = cosTheta > 0.f ? line + OMAT : line + IMAT ;
-    const int su_line = cosTheta > 0.f ? line + ISUR : line + OSUR ;
-
+    const int line = boundary * _BOUNDARY_NUM_MATSUR; // now that are not signing boundary use 0-based
+
+    int       m1_line = cosTheta > 0.f ? line + IMAT : line + OMAT;
+    const int m2_line = cosTheta > 0.f ? line + OMAT : line + IMAT;
+    const int su_line = cosTheta > 0.f ? line + ISUR : line + OSUR;
+
+    // A boundary normally describes a daughter volume relative to its parent.
+    // When two daughter volumes touch, however, a photon can move directly from
+    // one sibling to the other. The boundary selected by OptiX may then name the
+    // common parent instead of the material the photon is actually leaving.
+    //
+    // qsim carries the material-table row for the photon's current volume across
+    // interactions. Use that row as material1 when it is valid so absorption and
+    // scattering are sampled from the correct medium. Material2 and the surface
+    // still come from the selected boundary.
+    //
+    // UINT_MAX means that no current-material row is available. Line zero is a
+    // valid material row. Surface rows and rows beyond the uploaded boundary
+    // table must not be used for property or optical-buffer lookups.
+    const unsigned num_matline = boundary_meta->q0.u.y / _BOUNDARY_NUM_FLOAT4;
+    const unsigned matline_slot = carried_matline & 0x3u;
+    const bool     carried_matline_valid = carried_matline < num_matline && (matline_slot == unsigned(OMAT) || matline_slot == unsigned(IMAT));
+
+    if (carried_matline_valid)
+    {
+        m1_line = int(carried_matline);
+    }
 
-    s.material1 = boundary_lookup( wavelength, m1_line, 0);   // refractive_index, absorption_length, scattering_length, reemission_prob
-    s.m1group2  = boundary_lookup( wavelength, m1_line, 1);   // x: material1 group_velocity, y: material2 group_velocity, z/w unused
-    s.material2 = boundary_lookup( wavelength, m2_line, 0);   // refractive_index, (absorption_length, scattering_length, reemission_prob) only m2:refractive index actually used
-    s.surface   = boundary_lookup( wavelength, su_line, 0);   // detect,         , absorb            , (reflect_specular), reflect_diffuse     [they add to 1. so one not used]
+    s.material1 = boundary_lookup(wavelength, m1_line, 0); // refractive_index, absorption_length, scattering_length, reemission_prob
+    s.m1group2 = boundary_lookup(wavelength, m1_line, 1);  // x: material1 group_velocity, y: material2 group_velocity, z/w unused
+    s.material2 = boundary_lookup(wavelength, m2_line, 0); // refractive_index, (absorption_length, scattering_length, reemission_prob) only m2:refractive index actually used
+    s.surface = boundary_lookup(wavelength, su_line, 0);   // detect,         , absorb            , (reflect_specular), reflect_diffuse     [they add to 1. so one not used]
     s.set_material2_group_velocity(boundary_lookup(wavelength, m2_line, 1).x);
 
-    s.optical = optical[su_line].u ;   // 1-based-surface-index-0-meaning-boundary/type/finish/value  (type,finish,value not used currently)
+    s.optical = optical[su_line].u; // 1-based-surface-index-0-meaning-boundary/type/finish/value  (type,finish,value not used currently)
 
-    s.index.x = optical[m1_line].u.x ; // m1 index (1-based, see sstandard::make_optical)
-    s.index.y = optical[m2_line].u.x ; // m2 index (1-based, see sstandard::make_optical)
-    s.index.z = optical[su_line].u.x ; // su index (1-based, see sstandard::make_optical)
-    s.index.w = 0u ;                   // avoid undefined memory comparison issues
+    s.index.x = optical[m1_line].u.x; // m1 index (1-based, see sstandard::make_optical) -- reflects override if any
+    s.index.y = optical[m2_line].u.x; // m2 index (1-based, see sstandard::make_optical)
+    s.index.z = optical[su_line].u.x; // su index (1-based, see sstandard::make_optical)
+    s.index.w = 0u;                   // avoid undefined memory comparison issues
 
 #if !defined(PRODUCTION) && defined(DEBUG_PIDX)
-    if( idx == base_pidx )
+    if (idx == base_pidx)
     {
         printf("//qbnd.fill_state idx %7lld boundary %d line %d wavelength %10.4f m1_line %d m2_line %d su_line %d s.optical.x %d  \n", idx, boundary, line, wavelength, m1_line, m2_line, su_line, s.optical.x );
         printf("//qbnd.fill_state idx %7lld boundary %d [s.index.x-1](m1_index) %d [s.index.y-1](m2_index) %d [s.index.z-1](su_index) %d \n", idx, boundary, s.index.x-1u, s.index.y-1u, s.index.z-1u );