diff --git a/CSGOptiX/CSGOptiX.cc b/CSGOptiX/CSGOptiX.cc index 2b43e4e02..ae1a236f2 100644 --- a/CSGOptiX/CSGOptiX.cc +++ b/CSGOptiX/CSGOptiX.cc @@ -609,6 +609,59 @@ void CSGOptiX::initSimulate() 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 + + // Per-boundary face-bias gate for the sibling-coincident t-bias (see __intersection__is). + // bit=1 only where the boundary's own index contrast |n_m1 - n_m2| (max over wavelength) + // exceeds a TIR-risk threshold: high-contrast faces (quartz/air ~0.47) NEED the entering- + // sibling bias; low-contrast thin gaps (aerogel/air ~0.025) are harmed and gate OFF. + // Computed from the bnd table; no per-detector hardcode. Allocated once, reused per launch. + { + const SSim* ssim = foundry ? foundry->getSim() : nullptr; + const NP* bnd = ssim ? ssim->get_bnd() : nullptr; + if (bnd != nullptr && bnd->shape.size() == 5) + { + const unsigned num_bnd = bnd->shape[0]; + const unsigned num_wl = bnd->shape[3]; + const double* bv = bnd->cvalues(); // stree/standard bnd is bits(num_bnd, 0u); + for (unsigned b = 0; b < num_bnd; b++) + { + double dnmax = 0.; + for (unsigned w = 0; w < num_wl; w++) + { + const double n_om = bv[(((b * 4u + LMAT_OMAT) * 2u + 0u) * num_wl + w) * 4u + 0u]; // RINDEX of omat + const double n_im = bv[(((b * 4u + LMAT_IMAT) * 2u + 0u) * num_wl + w) * 4u + 0u]; // RINDEX of imat + double dn = n_om - n_im; + if (dn < 0.) + dn = -dn; + if (dn > dnmax) + dnmax = dn; + } + bits[b] = (dnmax > DN_THRESHOLD) ? 1u : 0u; + LOG(LEVEL) << "boundary_face_bias bnd " << b << " dnmax " << dnmax << " gate " << unsigned(bits[b]); + } + static unsigned char* s_face_bias = nullptr; + static unsigned s_face_bias_num = 0u; + if (s_face_bias == nullptr) + { + CUDA_CHECK(cudaMalloc((void**)&s_face_bias, num_bnd)); + s_face_bias_num = num_bnd; + } + if (s_face_bias_num == num_bnd) + { + CUDA_CHECK(cudaMemcpy(s_face_bias, bits.data(), num_bnd, cudaMemcpyHostToDevice)); + params->boundary_face_bias = s_face_bias; + } + LOG(LEVEL) << "boundary_face_bias gate computed, num_bnd " << num_bnd; + } + else + { + LOG(LEVEL) << "boundary_face_bias: no bnd table (rank != 5) -> gate left null (t-bias always-on)"; + } + } + params->tmax = 1000000.f ; params->max_time = SEventConfig::MaxTime() ; diff --git a/CSGOptiX/CSGOptiX7.cu b/CSGOptiX/CSGOptiX7.cu index 6fd39efd8..d462a0296 100644 --- a/CSGOptiX/CSGOptiX7.cu +++ b/CSGOptiX/CSGOptiX7.cu @@ -439,6 +439,14 @@ static __forceinline__ __device__ void simulate( const uint3& launch_idx, const #endif sim->generate_photon(ctx.p, rng, gs, photon_idx, genstep_idx ); + // F7: seed carried matline from the genstep's source-medium bnd row + // (gs.q0.u.z == scerenkov/sscint matline). Updated through + // propagate_at_boundary. Guard the -1 sentinel (0xFFFFFFFF) that + // SEvt::setGenstep emits for bad_ck/unmapped mtindex. + { + const unsigned gm = gs.q0.u.z; + ctx.current_matline = (gm == 0xFFFFFFFFu) ? 0u : gm; + } int command = START ; int bounce = 0 ; @@ -906,8 +914,19 @@ 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 ) ; + // Sibling-touching coincident-face fix: at a face shared by two sibling primitives OptiX + // sorts the two coplanar intersections inconsistently; the exiting hit (cosI>0) carries the + // parent medium as m2, wrong for a sibling crossing. Bias the t REPORTED to OptiX's closest- + // hit sort by +10 um on the exit side so the entering sibling wins; the real unbiased isect.w + // still goes to the PRD/attributes so position advancement is unperturbed. Gated per boundary + // by index contrast (params.boundary_face_bias from CSGOptiX::initSimulate): high-contrast + // faces (quartz/air ~0.47) need it; low-contrast thin gaps gate off; nullptr => always-on. + const float cosI = dot(ray_direction, make_float3(isect.x, isect.y, isect.z)); + const bool face_bias_on = (params.boundary_face_bias == nullptr) || (params.boundary_face_bias[boundary] != 0u); + const float t_report = (cosI > 0.f && face_bias_on) ? isect.w + 1.e-2f : isect.w; + #ifdef WITH_PRD - if(optixReportIntersection( isect.w, hitKind)) + if (optixReportIntersection(t_report, hitKind)) { quad2* prd = SOPTIX_getPRD(); // access prd addr from RG program prd->q0.f = isect ; // .w:distance and .xyz:normal which starts as the local frame one @@ -923,7 +942,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 diff --git a/CSGOptiX/Params.cc b/CSGOptiX/Params.cc index cc3a8b7c5..b703c27e4 100644 --- a/CSGOptiX/Params.cc +++ b/CSGOptiX/Params.cc @@ -126,9 +126,7 @@ std::string Params::detail() const 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), @@ -147,6 +145,7 @@ Params::Params(int raygenmode_, unsigned width, unsigned height, unsigned depth) origin_y(0), tmin(0.f), tmin0(0.f), + boundary_face_bias(nullptr), PropagateEpsilon0Mask(0u), tmax(0.f), vizmask(0xff), diff --git a/CSGOptiX/Params.h b/CSGOptiX/Params.h index 54a93faea..b99587e3f 100644 --- a/CSGOptiX/Params.h +++ b/CSGOptiX/Params.h @@ -64,6 +64,7 @@ struct Params float tmin ; float tmin0 ; + const unsigned char* boundary_face_bias; // per-boundary gate (1=apply sibling-coincident t-bias) keyed on |n_m1-n_m2|>0.1; nullptr => always-on unsigned PropagateEpsilon0Mask ; // default from SEventConfig TO,CK,SI,SC,RE float PropagateRefineDistance ; bool PropagateRefine ; diff --git a/qudarap/qbnd.h b/qudarap/qbnd.h index 57ff4a02c..9446a9585 100644 --- a/qudarap/qbnd.h +++ b/qudarap/qbnd.h @@ -38,7 +38,7 @@ struct qbnd #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 ); + 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 = 0u); #endif }; @@ -183,14 +183,30 @@ 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 ; + 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 ; + // F7 sibling-pair material override. + // GBndLib does not enumerate boundary entries for CSG sibling-pair faces + // (two daughters of one parent sharing a face), so OptiX may pick a boundary + // whose m1/m2 do not reference the photon's actual medium -> wrong material1. + // qsim::propagate passes the photon's carried matline (seeded from the + // genstep, updated each propagate_at_boundary); when it is valid we use it + // directly as m1 so absorption/scattering sample the true current medium. + // (material2 is left as the boundary's other side.) + // Guards: 0 = "no carry" (legacy callers); 0xFFFFFFFF = the lookup_mtline=-1 + // sentinel from SEvt::setGenstep, not a valid table row; and the carried + // slot must be a material row (OMAT/IMAT), not a surface row (OSUR/ISUR). + const unsigned matline_slot = carried_matline & 0x3u; + if (carried_matline != 0u && carried_matline != 0xFFFFFFFFu && (matline_slot == unsigned(OMAT) || matline_slot == unsigned(IMAT))) + { + 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 @@ -200,7 +216,7 @@ inline QBND_METHOD void qbnd::fill_state(sstate& s, unsigned boundary, float wav 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.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 diff --git a/qudarap/qsim.h b/qudarap/qsim.h index a4b7143d6..4968b7967 100644 --- a/qudarap/qsim.h +++ b/qudarap/qsim.h @@ -1263,6 +1263,19 @@ inline QSIM_METHOD int qsim::propagate_at_boundary(unsigned& flag, RNG& rng, sct ctx.current_material_index = reflect ? s.index.x : s.index.y; ctx.current_group_velocity = reflect ? s.material1_group_velocity() : s.material2_group_velocity(); + // F7 SIBLING-PAIR carry: update current_matline so propagate() at the NEXT + // boundary knows the photon's actual medium even when that boundary doesn't + // reference it (the sibling-pair case). _c1 = -dot(p.mom,normal) here, so + // _c1 < 0 maps to cosTheta > 0 in qbnd::fill_state (m1=IMAT, m2=OMAT). + { + const unsigned bnd_idx = ctx.prd->boundary(); + const unsigned imat_line = bnd_idx * _BOUNDARY_NUM_MATSUR + IMAT; + const unsigned omat_line = bnd_idx * _BOUNDARY_NUM_MATSUR + OMAT; + const unsigned this_m1_line = (_c1 < 0.f) ? imat_line : omat_line; + const unsigned this_m2_line = (_c1 < 0.f) ? omat_line : imat_line; + ctx.current_matline = reflect ? this_m1_line : this_m2_line; + } + #if !defined(PRODUCTION) && defined(DEBUG_TAG) if( flag == BOUNDARY_REFLECT ) { @@ -2153,7 +2166,7 @@ inline QSIM_METHOD int qsim::propagate(const int bounce, RNG& rng, sctx& ctx ) // copy geometry info into the sphoton struct ctx.p.set_prd(boundary, identity, cosTheta, iindex ); // HMM: lposcost not passed along - bnd->fill_state(ctx.s, boundary, ctx.p.wavelength, cosTheta, ctx.pidx, base->pidx ); + bnd->fill_state(ctx.s, boundary, ctx.p.wavelength, cosTheta, ctx.pidx, base->pidx, ctx.current_matline); // F7: pass carried matline for sibling-pair override if (ctx.current_material_index == 0u) { diff --git a/sysrap/sctx.h b/sysrap/sctx.h index 82e1dd389..f935577ab 100644 --- a/sysrap/sctx.h +++ b/sysrap/sctx.h @@ -79,6 +79,7 @@ struct sctx sstate s ; float current_group_velocity = 0.f; unsigned current_material_index = 0u; + unsigned current_matline = 0u; // F7 sibling-pair: photon's authoritative current-medium matline, seeded at birth, updated on transmit/reflect #ifndef PRODUCTION srec rec ;