// Full isotropic + uniaxial-birefringent + biaxial-birefringent spectral estimator, a
// direct translation of `optics::raytracer::trace_spectral_ray`. Driven by
// `renderer::gpu::estimator_check` (statistical image comparison, energy-conservation
// furnace anchor, spectral-space debug comparison) and `renderer::gpu::transport_check`
// (per-function ULP budgets, exercised via `shaders/transport_functions.wgsl`).
//
// "Mode A"/"mode B" generalize uniaxial ordinary/extraordinary to biaxial crystals with
// no single optic axis: mode A is the faster (lower-index) root of
// `biaxial_wave_indices`, mode B the slower; `is_extraordinary` selects between them and
// is reused unchanged from the uniaxial path. An air->crystal entry resolves both
// biaxial modes' wave-normal directions via `biaxial_resolve_entry_mode`'s two-iteration
// fixed point, evaluated once from the hero channel and shared by every companion
// channel. Both biaxial modes walk off (unlike uniaxial's ordinary ray, which never
// does). Whether this port is actually TRUSTED for a real render is governed entirely
// by `optics::materials::GemMaterial::gpu_supported` -- callers must consult that
// predicate, not infer support from this shader existing; see its own doc comment for
// current status.
//
// For an isotropic material, `is_anisotropic`/`is_biaxial` are always false by
// construction, which collapses this kernel back to bit-for-bit isotropic-only
// behaviour. The pleochroic Beer-Lambert absorption path is still ported faithfully
// (`electric_field_direction`, eigen-polarizations, the quadratic form) rather than
// shortcut to a bare `spectral_absorption(lambda)`, so even the isotropic limit
// reproduces the CPU's actual rounding, not merely its mathematically-equal result.
//
// Uniaxial interfaces route through the exact closed-form `entry_solve_pair`/
// `internal_solve` machinery (`transport_physics.wgsl`, mirroring
// `optics::raytracer::uniaxial_fresnel`), except at the exact degenerate
// wave-normal-parallel-to-optic-axis limit (`uniaxial_nondegenerate` false), where both
// eigenmodes truly collapse to a single isotropic response and the simpler
// scalar-Fresnel-per-mode path is exact rather than an approximation. Chromatic
// termination (a companion channel's own refracted/walk-off direction failing to match
// the shared hero-driven direction within `DIRECTION_MATCH_COS_TOL`) is reachable via
// either a dispersive Snell-angle or walk-off mismatch; the direction compared against
// must always be the stored `final_refr_dir`/`final_dir_k`, never recomputed a second
// time (a few-ULP-different recomputation of the same direction would fail its own
// match against itself). Verified against the CPU functions by
// `renderer::gpu::transport_check::p2_uniaxial_fresnel` and a Tier 3 image comparison
// on Zircon/Tourmaline/Quartz/Rutile.
//
// The pleochroic absorption path additionally consults the material's third
// (`beta_ray`) band set, when present, via `pleochroic_channel_alpha_biaxial`'s
// three-independent-coefficient quadratic form whenever `is_biaxial`.
//
// Every stochastic decision (Fresnel reflect/transmit, Russian roulette) uses its own
// locally computed probability for both the branch comparison and the compensating
// division, so float divergence between threads never affects correctness. Precision
// is `f32` throughout; `f32::mul_add` mirrors WGSL `fma`, and any power path uses a
// multiplication chain rather than `pow()`.
//
// `transport_main` always writes all four output buffers (final XYZ, and the
// pre-integration per-channel radiance/lambda/path_pdf debug arrays) rather than using
// two entry points with different binding subsets: WGSL's per-entry-point bind-group
// inference is based on static reachability, not runtime branching, so a `write_debug`
// parameter would still force both entry points to bind everything anyway.
// `estimator_check`'s large statistical dispatches simply allocate (never read back)
// the debug buffers.
//
// Per-thread state (Stokes vectors, path_pdf, lambdas, ray origin/dir, loop scalars)
// lives entirely in this function's local variables -- no cross-thread communication,
// no atomics, one thread per (pixel, sample) tuple writing only its own output slot.
// This is what makes two dispatches against identical input byte-identical
// (`estimator_check::run_determinism`): GPU scheduling order can never matter when no
// thread ever reads another thread's output.
//
// Exit-event spectral splitting: at a crystal->air exit, a still-alive companion
// channel `k` whose own Snell direction diverges from the hero's gets its own
// transmission and one bounded environment probe (`try_split_exit_channel`) instead of
// being zeroed outright, accumulating into a per-channel `split_radiance` that commits
// into `radiance` once, only if the shared path escapes. A per-channel `compat`
// bitmask tracks which channels remain in the same MIS family; an interior mismatch
// narrows it via `narrow_compat` (exits never narrow a family). Both Russian-roulette
// sites rescale `split_radiance` by the same `1/q` as `stokes` on survival. The final
// XYZ integration normalizes each channel over its own family via
// `optics::raytracer::color::integrate_channels_to_xyz_families`'s per-channel weight
// (`weight_k = N * path_pdf[hero] / sum_{j in compat[k]} path_pdf[j]`) rather than one
// shared weight. See each call site below for the exact CPU function it mirrors.
// Verified by `renderer::gpu::transport_check::p6_exit_splitting` and a Tier 3 image
// comparison isolating refractive (non-TIR-only) paths for a strongly dispersive
// uniaxial material.
//
// Kernel specialisation: `MATERIAL_CLASS` is a pipeline-overridable constant (WGSL
// `override`, resolved at pipeline-creation time via
// `wgpu::PipelineCompilationOptions::constants`) that lets `GpuFrameRenderer` compile
// one specialised pipeline per material class plus the generic pipeline every self-test
// uses unmodified (default 0u == generic, so nothing dispatching without the override
// is affected). It feeds only the `is_anisotropic`/`is_biaxial` derivations: forcing
// either false at compile time makes every per-ray array gated on it (biaxial axes,
// per-channel biaxial/uniaxial indices, wave directions, etc.) unreachable, so
// naga/the driver's dead-code elimination can drop the writes, the arrays, and the
// register space they would otherwise hold live across the whole bounce loop --
// isotropic and uniaxial pipelines never carry state they don't need.
// MATERIAL_CLASS_GENERIC (and MATERIAL_CLASS_BIAXIAL for `is_biaxial`) passes the
// runtime buffer value through unchanged, so the generic pipeline's behaviour is
// identical to before this feature existed.
override MATERIAL_CLASS: u32 = 0u;
const MATERIAL_CLASS_GENERIC: u32 = 0u;
const MATERIAL_CLASS_ISOTROPIC: u32 = 1u;
const MATERIAL_CLASS_UNIAXIAL: u32 = 2u;
const MATERIAL_CLASS_BIAXIAL: u32 = 3u;
const NUM_CHANNELS: u32 = 8u;
const SPECTRUM_MIN: f32 = 380.0;
const SPECTRUM_SPAN: f32 = 400.0;
const NORM_FACTOR: f32 = (400.0 / 8.0) / 106.856;
const DIRECTION_MATCH_COS_TOL: f32 = 1.0 - 1e-6;
const RR_FLOOR: f32 = 0.05;
// Reflect-vs-transmit SELECTION probability's own clamp bounds -- distinct from
// `R_UNPOL_MIN`/`R_UNPOL_MAX` (`transport_physics.wgsl`), which only ever scale
// `path_pdf`, never divide `stokes` directly. This megakernel's r_unpol entry decision
// is the one site that both drives a stochastic branch AND divides `stokes` by it
// (`1/r_unpol`, `1/(1-r_unpol)`), so it alone gets the tighter [0.02, 0.98] range --
// see `apply_partial_fresnel_bounce`'s doc comment on the CPU side.
const R_UNPOL_SELECT_MIN: f32 = 0.02;
const R_UNPOL_SELECT_MAX: f32 = 0.98;
// RAY_EPS, R_UNPOL_MIN, R_UNPOL_MAX, FRESNEL_BRANCH_STREAM, RUSSIAN_ROULETTE_STREAM,
// BIREFRINGENT_SPLIT_STREAM, MODE_COUPLING_STREAM, FROSTED_DIR_U_STREAM,
// FROSTED_DIR_V_STREAM, and hash_u32 live in `transport_physics.wgsl` so
// `apply_frosted_bounce` (shared with `transport_functions.wgsl`) has them in scope
// without a duplicate copy.
// optics::raytracer::{BRADFORD_XYZ_TO_LMS, BRADFORD_LMS_TO_XYZ, apply_von_kries_white_balance}.
// `params.white_balance` is `compute_illuminant_white_balance`'s precomputed
// Bradford-LMS-space scale -- applying it correctly means transforming to this same
// Bradford LMS basis, scaling, and transforming back, not multiplying into XYZ
// directly. Local to this file since `transport_functions.wgsl`'s Tier 2 kernels never
// apply a white balance -- only `compute_illuminant_white_balance` itself is tested,
// via `shaders/environment.wgsl`'s own separately defined copy of these constants.
const BRADFORD_XYZ_TO_LMS = mat3x3<f32>(
vec3<f32>(0.8951, -0.7502, 0.0389),
vec3<f32>(0.2664, 1.7135, -0.0685),
vec3<f32>(-0.1614, 0.0367, 1.0296),
);
const BRADFORD_LMS_TO_XYZ = mat3x3<f32>(
vec3<f32>(0.986993, 0.432305, -0.008529),
vec3<f32>(-0.147054, 0.518360, 0.040043),
vec3<f32>(0.159963, 0.049291, 0.968487),
);
fn apply_von_kries_white_balance(xyz: vec3<f32>, lms_scale: vec3<f32>) -> vec3<f32> {
let lms = BRADFORD_XYZ_TO_LMS * xyz;
return BRADFORD_LMS_TO_XYZ * (lms * lms_scale);
}
// ---------------------------------------------------------------------------------
// Struct layouts -- must match `renderer::buffers` field-for-field (see that module's
// doc comment on why a hand-derived offset is never trusted without the echo test).
// ---------------------------------------------------------------------------------
struct GpuCameraParams {
origin: vec3<f32>,
fov_tan: f32,
forward: vec3<f32>,
width: f32,
right: vec3<f32>,
height: f32,
up: vec3<f32>,
num_samples: u32,
}
struct GpuTransportParams {
num_pixels: u32,
max_bounces: u32,
sample_offset: u32,
env_mode: u32,
l0: f32,
studio_temp_k: f32,
studio_spot_mult: f32,
studio_exposure: f32,
studio_light_yaw: f32,
studio_light_pitch: f32,
pixel_offset: u32,
// Reused pad field: gates `transport_main`'s three per-channel debug writes.
write_debug_buffers: u32,
white_balance: vec3<f32>,
// Reused pad field: selects the tabulated CIE D65 measured spectrum over
// `blackbody_spectrum`.
studio_use_d65: u32,
}
struct DispersionParams {
model_type: u32,
param_a: vec4<f32>,
param_b: vec4<f32>,
param_c: vec4<f32>,
c_axis_and_birefringence: vec4<f32>,
is_anisotropic: u32,
biaxial_delta_beta_alpha: f32,
has_biaxial_delta: u32,
}
struct GpuGemMaterial {
dispersion: DispersionParams,
crystal_system: u32,
optical_character: u32,
is_pleochroic: u32,
o_ray_band_count: u32,
e_ray_band_count: u32,
o_ray_bands: array<AbsorptionBand, 8>,
e_ray_bands: array<AbsorptionBand, 8>,
scattering_sigma_s: f32,
scattering_g: f32,
edge_rounding_radius: f32,
// Appended fields below are always added at the end, never inserted earlier, so no
// existing field's offset ever shifts -- see `renderer::buffers::GpuGemMaterial`'s
// doc comment.
has_beta_ray: u32,
beta_ray_band_count: u32,
beta_ray_bands: array<AbsorptionBand, 8>,
absorption_path_scale: f32,
has_extraordinary_dispersion: u32,
extraordinary_model_type: u32,
extraordinary_param_a: vec4<f32>,
extraordinary_param_b: vec4<f32>,
}
struct FacetPlane {
normal: vec3<f32>,
d: f32,
}
@group(0) @binding(0) var<uniform> camera: GpuCameraParams;
@group(0) @binding(1) var<uniform> params: GpuTransportParams;
@group(0) @binding(2) var<storage, read> material: GpuGemMaterial;
@group(0) @binding(3) var<storage, read> planes: array<FacetPlane>;
@group(0) @binding(4) var<storage, read_write> out_xyz: array<f32>;
@group(0) @binding(5) var<storage, read_write> out_radiance: array<f32>;
@group(0) @binding(6) var<storage, read_write> out_lambdas: array<f32>;
@group(0) @binding(7) var<storage, read_write> out_path_pdf: array<f32>;
// `optics::raytracer::FacetFinish`, one entry per `planes[i]`, PARALLEL to `planes` (a
// separate binding, not a widened `FacetPlane` -- see `renderer::buffers::facet_finish`'s
// doc comment for why). `facet_finish::FROSTED` (1u) routes that facet's bounce through
// `apply_frosted_bounce` instead of the polished TIR/reflect/refract dispatch below; any
// other value (including an out-of-bounds index) is `facet_finish::POLISHED` (0u).
@group(0) @binding(8) var<storage, read> facet_finishes: array<u32>;
// One `compat[k]` narrowed-family bitmask per channel, per (pixel, sample) tuple --
// same `write_debug_buffers`-gated, self-test-only contract as `out_radiance`/
// `out_lambdas`/`out_path_pdf` above.
@group(0) @binding(9) var<storage, read_write> out_compat: array<u32>;
const FACET_FINISH_FROSTED: u32 = 1u;
// hash_u32 -- optics::raytracer::hash_u32, bit-exact. Defined in
// `transport_physics.wgsl` -- look there, not here.
//
// optics::raytracer::{low_discrepancy_base2, radical_inverse_base,
// cranley_patterson_rotate, PIXEL_JITTER_X_ROTATION_STREAM,
// PIXEL_JITTER_Y_ROTATION_STREAM, HERO_WAVELENGTH_ROTATION_STREAM}, bit-exact (see
// shaders/rng_equivalence.wgsl for the dedicated GPU/CPU RNG self-test).
//
// jx/jy/hero_rand use three different prime bases (2, 3, 5), not the same base rotated
// three ways -- measured: same-base pairing made variance worse for the
// highest-variance pixels.
const PIXEL_JITTER_X_ROTATION_STREAM: u32 = 0xA511E9B3u;
const PIXEL_JITTER_Y_ROTATION_STREAM: u32 = 0x63D81B23u;
const HERO_WAVELENGTH_ROTATION_STREAM: u32 = 0x1B873593u;
fn low_discrepancy_base2(n: u32) -> f32 {
return f32(reverseBits(n)) / 4294967296.0;
}
// optics::raytracer::radical_inverse_base -- general prime-base radical inverse (base 2
// uses the faster bit-reversal path above instead). Uses plain `+`/`*`/`/` (no `fma()`),
// matching the CPU side's non-fused `+=`/`/=` exactly.
fn radical_inverse_base(n_in: u32, base: u32) -> f32 {
var n = n_in;
var val: f32 = 0.0;
var inv_base: f32 = 1.0 / f32(base);
loop {
if (n == 0u) {
break;
}
let digit = n % base;
val = fma(f32(digit), inv_base, val);
inv_base = inv_base / f32(base);
n = n / base;
}
return val;
}
fn cranley_patterson_rotate(x: f32, offset: f32) -> f32 {
let sum = x + offset;
return sum - floor(sum);
}
// cie_1931_cmf -- color::cie1931::cie_1931_cmf, ported identically to
// shaders/environment.wgsl / shaders/furnace.wgsl.
fn cmf_lobe(x: f32, mu: f32, sigma_lo: f32, sigma_hi: f32) -> f32 {
var sigma: f32;
if (x < mu) {
sigma = sigma_lo;
} else {
sigma = sigma_hi;
}
let t = (x - mu) / sigma;
return exp(-0.5 * t * t);
}
fn cie_1931_cmf(l: f32) -> vec3<f32> {
let x_mid = fma(1.056, cmf_lobe(l, 599.8, 37.9, 31.0), 0.362 * cmf_lobe(l, 442.0, 16.0, 26.7));
let x = fma(0.065, -cmf_lobe(l, 501.1, 20.4, 26.2), x_mid);
let y = fma(0.286, cmf_lobe(l, 530.9, 16.3, 31.1), 0.821 * cmf_lobe(l, 568.8, 46.9, 40.5));
let z = fma(0.681, cmf_lobe(l, 459.0, 26.0, 13.8), 1.217 * cmf_lobe(l, 437.0, 11.8, 36.0));
return vec3<f32>(max(x, 0.0), max(y, 0.0), max(z, 0.0));
}
// optics::renderer::env_map_spectrum::rgb_to_spectral_radiance -- used for the
// direction-independent "uniform furnace" environment (env_mode == 0u): a grey
// `EnvironmentMap::uniform(w, h, [l0, l0, l0])` that is still wavelength-dependent, so
// this is deliberately not flattened to a bare `l0` return.
fn asymmetric_gaussian(x: f32, mu: f32, sigma_lo: f32, sigma_hi: f32) -> f32 {
var sigma: f32;
if (x < mu) {
sigma = sigma_lo;
} else {
sigma = sigma_hi;
}
let t = (x - mu) / sigma;
return exp(-0.5 * t * t);
}
fn rgb_to_spectral_radiance(r: f32, g: f32, b: f32, lambda_nm: f32) -> f32 {
let rc = max(r, 0.0);
let gc = max(g, 0.0);
let bc = max(b, 0.0);
return fma(
rc, asymmetric_gaussian(lambda_nm, 615.0, 45.0, 65.0),
fma(
gc, asymmetric_gaussian(lambda_nm, 545.0, 45.0, 45.0),
bc * asymmetric_gaussian(lambda_nm, 465.0, 40.0, 45.0),
),
);
}
// optics::raytracer::sample_studio_environment (+ optics::studio_rig::StudioRig) --
// ported identically to shaders/environment.wgsl.
fn powi_u(base: f32, exp: u32) -> f32 {
var result: f32 = 1.0;
var b: f32 = base;
var e: u32 = exp;
loop {
if (e == 0u) {
break;
}
if ((e & 1u) == 1u) {
result = result * b;
}
b = b * b;
e = e >> 1u;
}
return result;
}
fn blackbody_spectrum(lambda_nm: f32, temp_k: f32) -> f32 {
let t_k = max(temp_k, 1000.0);
let h_c_k: f32 = 14388000.0;
let exp_val = exp(min(h_c_k / (lambda_nm * t_k), 80.0));
let exp_560 = exp(min(h_c_k / (560.0 * t_k), 80.0));
let denom = max(exp_val - 1.0, 1e-6);
let denom_560 = max(exp_560 - 1.0, 1e-6);
let ratio = denom_560 / denom;
return clamp(powi_u(560.0 / lambda_nm, 5u) * ratio, 0.01, 20.0);
}
// optics::raytracer::environment::{CIE_D65_SPD_380_780_10NM, d65_relative_spectral_power}
// (CIE 15:2004) -- the "D65 Daylight" preset's real measured table, used instead of
// `blackbody_spectrum` above. A `const` array, not a uniform/storage buffer: 41 `f32`s
// is small enough to inline directly, needing no extra bind-group slot or upload.
const CIE_D65_SPD_380_780_10NM: array<f32, 41> = array<f32, 41>(
49.9755, 54.6482, 82.7549, 91.4860, 93.4318, 86.6823, 104.865, 117.008, 117.812, 114.861,
115.923, 108.811, 109.354, 107.802, 104.790, 107.689, 104.405, 104.046, 100.000, 96.3342,
95.7880, 88.6856, 90.0062, 89.5991, 87.6987, 83.2886, 83.6992, 80.0268, 80.2146, 82.2778,
78.2842, 69.7213, 71.6091, 74.3496, 61.6045, 69.8856, 75.0870, 63.5928, 46.4182, 66.8054,
63.3828,
);
// Ported op-for-op from `optics::raytracer::environment::d65_relative_spectral_power`:
// same 560nm-normalization (dividing by the table's `100.000` entry), same
// clamp-to-table-edge behaviour outside 380-780nm.
fn d65_relative_spectral_power(lambda_nm: f32) -> f32 {
let start_nm: f32 = 380.0;
let step_nm: f32 = 10.0;
let last_index: u32 = 40u; // CIE_D65_SPD_380_780_10NM.len() - 1
let clamped = clamp(lambda_nm, start_nm, fma(step_nm, f32(last_index), start_nm));
let position = (clamped - start_nm) / step_nm;
let index0 = min(u32(floor(position)), last_index - 1u);
let index1 = index0 + 1u;
let frac = position - f32(index0);
let v0 = CIE_D65_SPD_380_780_10NM[index0];
let v1 = CIE_D65_SPD_380_780_10NM[index1];
return fma(frac, v1 - v0, v0) / 100.0;
}
const RING_LIGHT_COUNT: u32 = 16u;
fn studio_rig_key_dir(light_yaw: f32, light_pitch: f32) -> vec3<f32> {
let cos_lp = cos(light_pitch);
let sin_lp = sin(light_pitch);
let cos_ly = cos(light_yaw);
let sin_ly = sin(light_yaw);
return normalize(vec3<f32>(cos_lp * sin_ly, sin_lp, cos_lp * cos_ly));
}
fn studio_rig_fill_dir(light_yaw: f32, light_pitch: f32) -> vec3<f32> {
let fill_yaw = fma(PI, 0.78, light_yaw);
let fill_pitch = clamp(light_pitch * 0.65, 0.15, 1.2);
return normalize(vec3<f32>(cos(fill_pitch) * sin(fill_yaw), sin(fill_pitch), cos(fill_pitch) * cos(fill_yaw)));
}
fn studio_rig_ring_dir(i: u32, light_yaw: f32, sin_lp: f32) -> vec3<f32> {
let angle = fma(f32(i), PI * 2.0 / f32(RING_LIGHT_COUNT), light_yaw);
return normalize(vec3<f32>(cos(angle) * 0.75, sin_lp * 0.8, sin(angle) * 0.75));
}
// `key_dir`/`fill_dir`/`sin_lp` (the `StudioRig`-equivalent quantities) are constant
// across an entire ray, so the caller (`transport_main`'s miss branch) computes them
// once before its `NUM_CHANNELS` loop and passes them in, rather than this function
// recomputing them on every per-channel call -- mirrors
// `optics::raytracer::accumulate_miss_radiance` building them once per ray.
fn sample_studio_environment_with_rig(
dir_in: vec3<f32>,
lambda_nm: f32,
key_dir: vec3<f32>,
fill_dir: vec3<f32>,
sin_lp: f32,
) -> f32 {
let d = normalize(dir_in);
// Mirrors `sample_studio_environment_with_rig`'s `LightingPreset::Daylight` branch.
var spec_power: f32;
if (params.studio_use_d65 != 0u) {
spec_power = d65_relative_spectral_power(lambda_nm);
} else {
spec_power = blackbody_spectrum(lambda_nm, params.studio_temp_k);
}
let bg_val = max(fma(0.012, fma(d.y, 0.5, 0.5), 0.015), 0.005) * params.studio_exposure;
var radiance = bg_val * spec_power;
let key_dot = max(dot(d, key_dir), 0.0);
if (key_dot > 0.0) {
let softbox = powi_u(key_dot, 28u) * 12.0 * params.studio_spot_mult * params.studio_exposure;
radiance = fma(softbox, spec_power, radiance);
}
let fill_dot = max(dot(d, fill_dir), 0.0);
if (fill_dot > 0.0) {
let fill = powi_u(fill_dot, 18u) * 4.5 * params.studio_exposure;
radiance = fma(fill, spec_power, radiance);
}
for (var i: u32 = 0u; i < RING_LIGHT_COUNT; i = i + 1u) {
let ring_dir = studio_rig_ring_dir(i, params.studio_light_yaw, sin_lp);
let ring_dot = max(dot(d, ring_dir), 0.0);
if (ring_dot > 0.96) {
let spark = (ring_dot - 0.96) / 0.04;
let intensity = powi_u(spark, 6u) * 22.0 * params.studio_spot_mult * params.studio_exposure;
radiance = fma(intensity, spec_power, radiance);
}
}
return radiance;
}
fn sample_environment_with_rig(
dir: vec3<f32>,
lambda_nm: f32,
key_dir: vec3<f32>,
fill_dir: vec3<f32>,
sin_lp: f32,
) -> f32 {
if (params.env_mode == 0u) {
return rgb_to_spectral_radiance(params.l0, params.l0, params.l0, lambda_nm);
}
return sample_studio_environment_with_rig(dir, lambda_nm, key_dir, fill_dir, sin_lp);
}
// `dispersion_evaluate`, `spectral_absorption`, the four `mueller_*` Mueller-matrix
// constructors, `tir_phase_delta`, `normalize_or_zero`, `signed_frame_rotation_psi`,
// `degree_of_polarization`, `polarization_azimuth`, `arbitrary_perpendicular`,
// `electric_field_direction`, `stable_orthonormal_basis_t`,
// `ordinary_eigen_polarization`, `extraordinary_eigen_polarization`, `quadratic_form`,
// and `pleochroic_channel_alpha` all live in `shaders/transport_physics.wgsl`, the
// shared source `build.rs` concatenates ahead of this file. Look there, not here.
// optics::raytracer::intersect_polyhedron
struct HitInfo {
hit: bool,
t: f32,
normal: vec3<f32>,
// optics::raytracer::HitRecord::facet_idx -- lets the frosted-finish lookup below
// index `facet_finishes`. Mirrors `near_facet`/`far_facet`'s `.unwrap_or(0)`
// fallback: `0u` whenever a hit is never reported either.
facet_idx: u32,
}
fn intersect_ray(origin: vec3<f32>, dir: vec3<f32>) -> HitInfo {
var t_near: f32 = -1e30;
var t_far: f32 = 1e30;
var near_normal = vec3<f32>(0.0, 0.0, 0.0);
var far_normal = vec3<f32>(0.0, 0.0, 0.0);
var near_idx: u32 = 0u;
var far_idx: u32 = 0u;
var result: HitInfo;
let num_planes = arrayLength(&planes);
for (var i: u32 = 0u; i < num_planes; i = i + 1u) {
let p = planes[i];
let n = p.normal;
let denom = dot(n, dir);
let side = p.d + dot(n, origin);
let numer = -side;
if (abs(denom) > 1e-7) {
let t = numer / denom;
if (denom < 0.0) {
if (t > t_near) {
t_near = t;
near_normal = n;
near_idx = i;
}
} else if (t < t_far) {
t_far = t;
far_normal = n;
far_idx = i;
}
} else if (side > 0.0) {
// Ray (near-)parallel to this plane, origin already outside its half-space
// -- the polyhedron intersection is empty for this ray.
result.hit = false;
result.t = 0.0;
result.normal = vec3<f32>(0.0, 0.0, 0.0);
result.facet_idx = 0u;
return result;
}
}
if (t_near > t_far) {
result.hit = false;
result.t = 0.0;
result.normal = vec3<f32>(0.0, 0.0, 0.0);
result.facet_idx = 0u;
} else if (t_near > 1e-4) {
result.hit = true;
result.t = t_near;
result.normal = near_normal;
result.facet_idx = near_idx;
} else if (t_far > 1e-4) {
result.hit = true;
result.t = t_far;
result.normal = far_normal;
result.facet_idx = far_idx;
} else {
result.hit = false;
result.t = 0.0;
result.normal = vec3<f32>(0.0, 0.0, 0.0);
result.facet_idx = 0u;
}
return result;
}
// Exit-event spectral splitting: optics::raytracer::refraction::try_split_exit_channel
// -- the one bounded fan-out primitive every exit-event mismatch site below calls.
// Reads the `planes` binding (via `intersect_ray`) and `params`/`material` (via
// `sample_environment_with_rig`), so -- like `intersect_ray`/`shading_normal_near_edge`
// above -- this stays megakernel-local rather than living in the shared
// `transport_physics.wgsl` prelude. `split_radiance` is threaded as a pointer to the
// caller's own local array, like `stokes`/`path_pdf`.
fn try_split_exit_channel(
split_radiance: ptr<function, array<f32, 8>>,
hit_point: vec3<f32>,
k: u32,
lambda_k: f32,
dir_k: vec3<f32>,
transmitted_intensity: f32,
key_dir: vec3<f32>,
fill_dir: vec3<f32>,
sin_lp: f32,
) {
let probe = intersect_ray(hit_point + dir_k * 1e-4, dir_k);
if (probe.hit) {
// Bounded re-entry: decline to trace further (a pure energy-loss truncation).
return;
}
let env_spectral = sample_environment_with_rig(dir_k, lambda_k, key_dir, fill_dir, sin_lp);
(*split_radiance)[k] = fma(max(transmitted_intensity, 0.0), env_spectral, (*split_radiance)[k]);
}
// optics::raytracer::shading_normal_near_edge -- reads the `planes` storage binding
// directly (like `intersect_ray` above), so stays megakernel-local rather than living
// in the shared `transport_physics.wgsl` prelude. `shaders/shading_normal.wgsl` has its
// own standalone copy with its own `planes` binding for
// `renderer::gpu::transport_check::run_shading_normal_near_edge`'s Tier 2 self-test;
// both are unmodified line-for-line translations of the same CPU function.
fn shading_normal_near_edge(hit_point: vec3<f32>, hit_facet_idx: u32, hit_normal: vec3<f32>, rounding_radius: f32) -> vec3<f32> {
if (rounding_radius <= 0.0) {
return hit_normal;
}
var nearest_dist: f32 = 1e30;
var nearest_normal = hit_normal;
let num_planes = arrayLength(&planes);
for (var i: u32 = 0u; i < num_planes; i = i + 1u) {
if (i == hit_facet_idx) {
continue;
}
let p = planes[i];
let dist = -(p.d + dot(p.normal, hit_point));
if (dist < nearest_dist) {
nearest_dist = dist;
nearest_normal = p.normal;
}
}
if (nearest_dist >= rounding_radius) {
return hit_normal;
}
let t = clamp(1.0 - nearest_dist / rounding_radius, 0.0, 1.0);
let smooth_t = t * t * fma(-2.0, t, 3.0);
let bisector = normalize_or_zero(hit_normal + nearest_normal);
return normalize_or_zero(hit_normal * (1.0 - smooth_t) + bisector * smooth_t);
}
// ---------------------------------------------------------------------------------
// optics::raytracer::Camera::generate_ray -- bit-exact per Phase 1's camera_check.
// ---------------------------------------------------------------------------------
struct RayGen {
origin: vec3<f32>,
dir: vec3<f32>,
}
fn generate_camera_ray(pixel: u32, jx: f32, jy: f32) -> RayGen {
let x = f32(pixel % u32(camera.width));
let y = f32(pixel / u32(camera.width));
let aspect = camera.width / camera.height;
let u = ((x + jx) / camera.width - 0.5) * 2.0 * aspect * camera.fov_tan;
let v = (0.5 - (y + jy) / camera.height) * 2.0 * camera.fov_tan;
let dir = normalize(camera.forward + camera.right * u + camera.up * v);
var result: RayGen;
result.origin = camera.origin;
result.dir = dir;
return result;
}
// optics::raytracer::apply_internal_mode_coupling -- stochastic o<->e re-coupling at an
// INTERNAL reflection inside an anisotropic crystal. This is a RELABELING, not a SPLIT:
// unlike the entry split (BIREFRINGENT_SPLIT_STREAM, where one incident ray genuinely
// becomes two physically distinct rays each carrying its own ~0.5 energy share already
// matched by its own 0.5 selection probability, needing no `1/p` scaling), this draw
// only re-rolls which eigenmode label governs the index used for the path's next
// bounce -- no new ray, no energy split, nothing for `stokes`/`path_pdf` to do. This
// draw's own selection probability is polarization-weighted
// (`entry_eigenmode_selection`) rather than a blanket 50/50, which is safe precisely
// because `stokes`/`path_pdf` are never touched here.
// optics::raytracer::refraction::entry_eigenmode_selection -- polarization-weighted
// probability that the uniaxial ORDINARY eigenmode is the physically correct label for
// the hero-driven path, plus that eigenmode's own doubled polarization azimuth
// (`cos(2*psi_o)`, `sin(2*psi_o)`) in the SAME Stokes reference frame
// `current_plane_normal` establishes (Malus's law energy split; a unit-vector
// double-angle identity in place of an atan2/half-angle round-trip). WGSL has no
// `Option`, so `valid == 0u` stands in for the Rust `None` case (negligible `i`,
// degenerate plane of incidence, or negligible linear polarization) -- callers must
// check it before trusting `p_o`/`cos_2psi_o`/`sin_2psi_o`.
struct EntryEigenmodeSelection {
valid: u32,
p_o: f32,
cos_2psi_o: f32,
sin_2psi_o: f32,
}
fn entry_eigenmode_selection(
c_axis: vec3<f32>,
current_plane_normal: vec3<f32>,
k_hat: vec3<f32>,
hero_i: f32,
hero_q: f32,
hero_u: f32,
) -> EntryEigenmodeSelection {
var result: EntryEigenmodeSelection;
result.valid = 0u;
result.p_o = 0.5;
result.cos_2psi_o = 0.0;
result.sin_2psi_o = 0.0;
if (hero_i <= 1e-7 || dot(current_plane_normal, current_plane_normal) <= 1e-6) {
return result;
}
if (fma(hero_q, hero_q, hero_u * hero_u) <= 1e-12) {
return result;
}
let o_hat = ordinary_eigen_polarization(k_hat, c_axis);
// `current_plane_normal` is already unit and exactly perpendicular to `k_hat` --
// no re-orthogonalization needed.
let s_hat = current_plane_normal;
let p_hat = cross(k_hat, s_hat);
let s_comp = dot(o_hat, s_hat);
let p_comp = dot(o_hat, p_hat);
let cos_2psi_o = fma(s_comp, s_comp, -(p_comp * p_comp));
let sin_2psi_o = 2.0 * s_comp * p_comp;
result.valid = 1u;
result.p_o = clamp(0.5 + 0.5 * fma(hero_q, cos_2psi_o, hero_u * sin_2psi_o) / hero_i, 0.0, 1.0);
result.cos_2psi_o = cos_2psi_o;
result.sin_2psi_o = sin_2psi_o;
return result;
}
// `internal_mode_coupling_draw`'s own selection probability is weighted by
// `entry_eigenmode_selection` (evaluated at this bounce's post-reflection wave normal
// `new_k` and Stokes state) instead of a blanket 50/50 -- safe because this draw never
// touches `stokes`/`path_pdf` either way. `is_biaxial` (no uniaxial ordinary eigenaxis
// to weight against) keeps the blanket 50/50 unconditionally, matching
// `apply_internal_mode_coupling`'s own biaxial exclusion. `has_exact_p_o`/`exact_p_o`
// mirror `apply_internal_mode_coupling`'s `exact_p_o: Option<f32>` parameter (WGSL has
// no `Option`, so a presence flag plus a value stands in): the same closed-form
// `internal_solve` Poynting-weighted `R_o/(R_o+R_e)` split the uniaxial TIR-forced and
// general-internal-reflect dispatch arms below already compute for their own
// reflected-energy accounting, reused here instead of the polarization-projection
// heuristic. `has_exact_p_o == false` (every FROSTED/REFLECT-arm call site for an
// isotropic or biaxial material) preserves the pre-existing heuristic exactly.
fn internal_mode_coupling_draw(
c_axis: vec3<f32>,
is_biaxial: bool,
current_plane_normal: vec3<f32>,
new_k: vec3<f32>,
hero_i: f32,
hero_q: f32,
hero_u: f32,
has_exact_p_o: bool,
exact_p_o: f32,
seed0: u32,
bounce: u32,
) -> bool {
var p_o = 0.5;
if (has_exact_p_o) {
p_o = exact_p_o;
} else if (!is_biaxial) {
let sel = entry_eigenmode_selection(c_axis, current_plane_normal, new_k, hero_i, hero_q, hero_u);
if (sel.valid != 0u) {
p_o = sel.p_o;
}
}
let split_rand = f32(hash_u32(seed0 ^ hash_u32(bounce ^ MODE_COUPLING_STREAM))) / 4294967295.0;
return split_rand >= p_o;
}
// optics::raytracer::trace_spectral_ray -- the isotropic-only estimator.
@compute @workgroup_size(64)
fn transport_main(@builtin(global_invocation_id) gid: vec3<u32>) {
let idx = gid.x;
let total = params.num_pixels * camera.num_samples;
if (idx >= total) {
return;
}
// GLOBAL pixel index: `params.pixel_offset` is 0 for a whole-frame dispatch, non-zero
// only for `renderer::gpu::frame`'s chunked dispatches. Added HERE, not to `idx`, so
// it reaches camera-ray generation and the per-pixel Cranley-Patterson rotations
// below (a function of the pixel's place in the frame), while output writes at the
// end of this function still index by the dispatch-local `idx`.
let pixel = idx / camera.num_samples + params.pixel_offset;
let local_sample = idx % camera.num_samples;
let sample_num = local_sample + params.sample_offset;
let seed0 = hash_u32((pixel * 0x9e3779b9u) ^ (sample_num * 0x85ebca6bu));
// Stratified pixel jitter and hero wavelength (optics::raytracer::
// {low_discrepancy_base2, cranley_patterson_rotate}), not an unstratified
// hash-uniform. `seed0` above still seeds every per-bounce draw below (Fresnel
// branch, Russian roulette, birefringent split); only jx/jy/hero_rand come from
// this construction.
let rot_jx = low_discrepancy_base2(hash_u32(pixel ^ PIXEL_JITTER_X_ROTATION_STREAM));
let rot_jy = low_discrepancy_base2(hash_u32(pixel ^ PIXEL_JITTER_Y_ROTATION_STREAM));
let rot_hero = low_discrepancy_base2(hash_u32(pixel ^ HERO_WAVELENGTH_ROTATION_STREAM));
let jx = cranley_patterson_rotate(low_discrepancy_base2(sample_num), rot_jx) - 0.5;
let jy = cranley_patterson_rotate(radical_inverse_base(sample_num, 3u), rot_jy) - 0.5;
let hero_rand = cranley_patterson_rotate(radical_inverse_base(sample_num, 5u), rot_hero);
let raygen = generate_camera_ray(pixel, jx, jy);
let channel_width = SPECTRUM_SPAN / f32(NUM_CHANNELS);
let lambda_hero = fma(hero_rand, SPECTRUM_SPAN, SPECTRUM_MIN);
var lambdas: array<f32, 8>;
for (var k: u32 = 0u; k < NUM_CHANNELS; k = k + 1u) {
let offset = fma(f32(k), channel_width, lambda_hero - SPECTRUM_MIN);
lambdas[k] = SPECTRUM_MIN + (offset % SPECTRUM_SPAN);
}
var stokes: array<vec4<f32>, 8>;
var radiance: array<f32, 8>;
var path_pdf: array<f32, 8>;
for (var k: u32 = 0u; k < NUM_CHANNELS; k = k + 1u) {
stokes[k] = vec4<f32>(1.0, 0.0, 0.0, 0.0);
radiance[k] = 0.0;
path_pdf[k] = 1.0;
}
let c_axis = material.dispersion.c_axis_and_birefringence.xyz;
let birefringence_delta = material.dispersion.c_axis_and_birefringence.w;
// ANDed with the pipeline-overridable MATERIAL_CLASS (see this file's header) -- a
// no-op for MATERIAL_CLASS_GENERIC, forced `false` for MATERIAL_CLASS_ISOTROPIC
// regardless of the material buffer's own flag.
let is_anisotropic = (material.dispersion.is_anisotropic != 0u)
&& (MATERIAL_CLASS != MATERIAL_CLASS_ISOTROPIC);
// Whether this material's anisotropy is genuinely biaxial (three distinct principal
// indices) rather than the uniaxial ordinary/extraordinary approximation --
// hoisted per-ray, constant across every bounce, like `c_axis` above. Never true
// for MATERIAL_CLASS_ISOTROPIC/UNIAXIAL's routed materials anyway (see
// `renderer::gpu::frame::classify_material`), so ANDing with MATERIAL_CLASS only
// removes dead-by-construction reachability.
let is_biaxial = (material.dispersion.has_biaxial_delta != 0u)
&& (MATERIAL_CLASS == MATERIAL_CLASS_GENERIC || MATERIAL_CLASS == MATERIAL_CLASS_BIAXIAL);
// The biaxial principal-axis frame (alpha, beta, gamma world directions) depends
// only on `c_axis`, so it is computed once per ray rather than every bounce.
let biax_axes = biaxial_axes_from_gamma(c_axis);
let biax_ax0 = biax_axes.ax0;
let biax_ax1 = biax_axes.ax1;
let biax_ax2 = biax_axes.ax2;
// The hero channel's base dispersion value and, from it, the hero indicatrix's
// three principal indices (`optics::materials::GemMaterial::biaxial_indicatrix`'s
// convention: `n_beta := dispersion.evaluate(lambda)`, `n_alpha := n_beta -
// biaxial_delta_beta_alpha`, `n_gamma := n_alpha + birefringence_delta`). Computed
// unconditionally: harmless when `!is_biaxial` since `biaxial_delta_beta_alpha` is
// `0.0` for every non-biaxial material, so `n_alpha_hero == n_beta_hero ==
// n_o_hero_seed`, and no `biaxial_*` function is called with these values unless
// `is_biaxial` guards it.
let n_o_hero_seed = dispersion_evaluate(material.dispersion.model_type, material.dispersion.param_a, material.dispersion.param_b, lambdas[0]);
let n_beta_hero = n_o_hero_seed;
let n_alpha_hero = n_beta_hero - material.dispersion.biaxial_delta_beta_alpha;
let n_gamma_hero = n_alpha_hero + birefringence_delta;
// `n_o_ch[k]` and the per-channel pleochroic absorption coefficients depend only on
// `lambdas[k]` (fixed for the whole ray) and the material's dispersion/band data,
// never on anything that varies per bounce -- so hoisting them to per-ray arrays,
// computed once, is bit-identical to the previous per-bounce recomputation.
// `alpha_beta_hoisted` is left zero unless `is_biaxial`, avoiding pointless
// band-summation work for the common non-biaxial case.
//
// Hoisting measured substantially faster than per-bounce recomputation
// (best-of-3, 800x600 Zircon/Tourmaline, 24 bounces). Kept.
var n_o_hoisted: array<f32, 8>;
var alpha_o_hoisted: array<f32, 8>;
var alpha_e_hoisted: array<f32, 8>;
var alpha_beta_hoisted: array<f32, 8>;
for (var k: u32 = 0u; k < NUM_CHANNELS; k = k + 1u) {
n_o_hoisted[k] = dispersion_evaluate(material.dispersion.model_type, material.dispersion.param_a, material.dispersion.param_b, lambdas[k]);
alpha_o_hoisted[k] = spectral_absorption(material.o_ray_bands, material.o_ray_band_count, lambdas[k]);
alpha_e_hoisted[k] = spectral_absorption(material.e_ray_bands, material.e_ray_band_count, lambdas[k]);
alpha_beta_hoisted[k] = 0.0;
}
if (is_biaxial) {
for (var k: u32 = 0u; k < NUM_CHANNELS; k = k + 1u) {
alpha_beta_hoisted[k] = spectral_absorption(material.beta_ray_bands, material.beta_ray_band_count, lambdas[k]);
}
}
var current_origin = raygen.origin;
var current_dir = raygen.dir;
// The wave normal `k`, tracked alongside `current_dir` (the Poynting/energy
// direction `S`). Starts equal to the ray's direction: outside the gem (air)
// `k == S` always (isotropic medium, no walk-off).
var current_k = raygen.dir;
var inside_gem = false;
// Which eigenmode the ray currently inside the crystal was stochastically assigned
// to at its most recent air->crystal entry. Meaningless while `!inside_gem`;
// carried across internal bounces exactly as the CPU carries it.
var is_extraordinary = false;
var prev_plane_normal = vec3<f32>(0.0, 0.0, 0.0);
var have_prev_plane_normal = false;
// Studio rig key/fill/ring directions, built once per ray: a pure function of
// `params` alone, needed both by the final miss branch and by
// `try_split_exit_channel`'s environment probe at every exit event. Unused when
// `params.env_mode == 0u` (HDR/flat environment).
let studio_key_dir = studio_rig_key_dir(params.studio_light_yaw, params.studio_light_pitch);
let studio_fill_dir = studio_rig_fill_dir(params.studio_light_yaw, params.studio_light_pitch);
let studio_sin_lp = sin(params.studio_light_pitch);
// Exit-event spectral splitting's per-trace state (see this file's header).
// `split_radiance` is a staging accumulator folded into `radiance` only if
// `path_escaped` ends up true; `compat` tracks which channels have matched at
// every interior dispersive event so far (bit `j` of `compat[k]` set while `k`/`j`
// remain compatible; all bits set before the first mismatch).
var split_radiance: array<f32, 8>;
var compat: array<u32, 8>;
for (var k: u32 = 0u; k < NUM_CHANNELS; k = k + 1u) {
split_radiance[k] = 0.0;
compat[k] = 0xFFu;
}
var path_escaped = false;
for (var bounce: u32 = 0u; bounce < params.max_bounces; bounce = bounce + 1u) {
let hit = intersect_ray(current_origin, current_dir);
if (!hit.hit) {
for (var k: u32 = 0u; k < NUM_CHANNELS; k = k + 1u) {
let env = sample_environment_with_rig(current_dir, lambdas[k], studio_key_dir, studio_fill_dir, studio_sin_lp);
// `max(stokes[k].x, 0.0)` clamps `I` to >= 0 before the environment
// lookup, matching `accumulate_miss_radiance`'s `StokesVector::intensity`
// on the CPU side -- negative `I` is unphysical on either side.
radiance[k] = fma(max(stokes[k].x, 0.0), env, radiance[k]);
}
// The only site that sets this: staged exit-split contributions commit
// only when the shared/hero path itself reaches this environment lookup.
path_escaped = true;
break;
}
// Mirrors optics::raytracer::trace_spectral_ray_inner's restructured bounce
// loop: attempt a Henyey-Greenstein scattering event somewhere along this
// segment before the plane-of-incidence rotation / facet processing below.
// Gated on `material.scattering_sigma_s > 0.0` -- every scene with it `<= 0.0`
// skips this block entirely, matching the CPU's default-off bit-identity
// guarantee.
if (inside_gem && material.scattering_sigma_s > 0.0) {
var s_axis = vec3<f32>(0.0, 0.0, 0.0);
if (have_prev_plane_normal) {
s_axis = prev_plane_normal;
}
// Mirrors the absorption block's `is_biaxial` branching below:
// `try_scatter_step` feeds the same channel_absorption_alphas the
// absorption block uses, so a biaxial material's scattering alphas must
// come from the biaxial indicatrix's own eigen_polarizations (at the hero
// wavelength), not the uniaxial approximation. Guarded on `is_biaxial`
// (false for every non-biaxial scene), so bit-identical otherwise.
// Eigen-polarizations and the propagation direction fed to
// pleochroic_channel_alpha(_biaxial) are properties of the WAVE NORMAL `k`,
// not the Poynting direction `S` -- current_k, not current_dir.
var eigen_a: vec3<f32>;
var eigen_b: vec3<f32>;
if (is_biaxial) {
let eig = biaxial_eigen_polarizations(n_alpha_hero, n_beta_hero, n_gamma_hero, biax_ax0, biax_ax1, biax_ax2, current_k);
eigen_a = eig.d_slow;
eigen_b = eig.d_fast;
} else {
eigen_a = ordinary_eigen_polarization(current_k, c_axis);
eigen_b = extraordinary_eigen_polarization(current_k, c_axis);
}
var alphas: array<f32, 8>;
for (var k: u32 = 0u; k < NUM_CHANNELS; k = k + 1u) {
// alpha_o/alpha_e/alpha_beta hoisted above (per-ray, not per-bounce).
if (is_biaxial && material.has_beta_ray != 0u) {
alphas[k] = pleochroic_channel_alpha_biaxial(
alpha_o_hoisted[k], alpha_beta_hoisted[k], alpha_e_hoisted[k], c_axis, s_axis, current_k, eigen_a, eigen_b, stokes[k],
);
} else {
alphas[k] = pleochroic_channel_alpha(
alpha_o_hoisted[k], alpha_e_hoisted[k], c_axis, s_axis, current_k, eigen_a, eigen_b, stokes[k],
);
}
}
let sc = maybe_scatter_or_extinguish(
alphas, material.scattering_sigma_s, material.scattering_g, current_dir, hit.t,
material.absorption_path_scale, seed0, bounce, &stokes, &path_pdf,
);
if (sc.scattered != 0u) {
current_origin = current_origin + sc.t_free * current_dir;
current_dir = sc.new_dir;
// A scattering event depolarizes, so `k` collapses to `S` going forward.
current_k = sc.new_dir;
// Scattered Stokes vectors are already depolarized, so the previous
// plane of incidence is no longer physically meaningful -- reset it
// like the pre-first-bounce state.
have_prev_plane_normal = false;
if (bounce > 4u) {
var max_intensity: f32 = 0.0;
for (var k: u32 = 0u; k < NUM_CHANNELS; k = k + 1u) {
max_intensity = max(max_intensity, max(stokes[k].x, 0.0));
}
let q = clamp(max_intensity, RR_FLOOR, 1.0);
let rr_rand = f32(hash_u32(seed0 ^ hash_u32(bounce ^ RUSSIAN_ROULETTE_STREAM))) / 4294967295.0;
if (rr_rand > q) {
break;
}
// `split_radiance` rides along on the same `1/q` survival rescale
// as `stokes`.
for (var k: u32 = 0u; k < NUM_CHANNELS; k = k + 1u) {
stokes[k] = stokes[k] * (1.0 / q);
split_radiance[k] = split_radiance[k] / q;
}
}
continue;
}
// No scatter event fired: the path survived to the facet boundary, and
// `maybe_scatter_or_extinguish` already applied this segment's full
// extinction weight (absorption AND the scattering removal probability)
// to `stokes`/`path_pdf`. Fall through to the facet-dispatch code below,
// but see the (skipped) absorption block further down for why it is not
// applied a second time.
}
let hit_point = current_origin + hit.t * current_dir;
// Task 2 (facet edge rounding): see shading_normal_near_edge's own doc comment.
var normal = shading_normal_near_edge(hit_point, hit.facet_idx, hit.normal, material.edge_rounding_radius);
// P2: `wave_dir_at_bounce` is the wave normal `k` (== current_k), used for
// EVERY index lookup, cos_i/sin_i, Snell/Fresnel evaluation, TIR decision, and
// the Stokes plane-of-incidence frame below -- see
// optics::raytracer::refraction's "wave normal vs Poynting direction" design
// note. The Poynting/energy direction `S` (== current_dir) is used directly
// (geometric origin-advance/intersection only, unaffected by P2) where still
// needed. `wave_dir_at_bounce == current_dir` trivially outside the crystal and
// for the uniaxial ordinary eigenmode, so every such case is bit-identical to
// the pre-P2 all-`current_dir` code.
let wave_dir_at_bounce = current_k;
// Plane-of-incidence frame rotation (Fix 1: signed psi via atan2).
let cpn_raw = cross(wave_dir_at_bounce, normal);
let cpn_len2 = dot(cpn_raw, cpn_raw);
var current_plane_normal: vec3<f32>;
if (cpn_len2 > 0.0) {
current_plane_normal = cpn_raw / sqrt(cpn_len2);
} else {
current_plane_normal = vec3<f32>(0.0, 0.0, 0.0);
}
if (have_prev_plane_normal && cpn_len2 > 1e-6 && dot(prev_plane_normal, prev_plane_normal) > 1e-6) {
let psi = signed_frame_rotation_psi(prev_plane_normal, current_plane_normal, wave_dir_at_bounce);
let rot = mueller_frame_rotation(psi);
for (var k: u32 = 0u; k < NUM_CHANNELS; k = k + 1u) {
stokes[k] = rot * stokes[k];
}
}
prev_plane_normal = current_plane_normal;
have_prev_plane_normal = true;
if (inside_gem) {
normal = -normal;
// Task 1: a scattering-active material's extinction for this segment was
// already applied above (in the `maybe_scatter_or_extinguish` no-scatter
// branch) -- applying the plain absorption loop here too would charge
// this segment's absorption TWICE. `material.scattering_sigma_s <= 0.0`
// is exactly the pre-Task-1 case (including every existing scene), where
// the block above never ran and this is the ONLY absorption
// application, matching the CPU's identical guard.
if (material.scattering_sigma_s <= 0.0) {
// Phase 4: optics::raytracer::channel_absorption_alphas -- for a
// biaxial material the two eigenmode directions come from the
// material's OWN indicatrix (`biaxial_eigen_polarizations`, evaluated
// at the hero wavelength -- the direction-only quantities this feeds
// don't vary per channel, exactly like the uniaxial ordinary/
// extraordinary pair they generalize), not the uniaxial
// ordinary/extraordinary approximation.
var eigen_a: vec3<f32>;
var eigen_b: vec3<f32>;
if (is_biaxial) {
let eig = biaxial_eigen_polarizations(n_alpha_hero, n_beta_hero, n_gamma_hero, biax_ax0, biax_ax1, biax_ax2, wave_dir_at_bounce);
eigen_a = eig.d_slow;
eigen_b = eig.d_fast;
} else {
eigen_a = ordinary_eigen_polarization(wave_dir_at_bounce, c_axis);
eigen_b = extraordinary_eigen_polarization(wave_dir_at_bounce, c_axis);
}
for (var k: u32 = 0u; k < NUM_CHANNELS; k = k + 1u) {
// R3-GPU: alpha_o/alpha_e/alpha_beta hoisted above (per-ray, not
// per-bounce) -- see that hoist's own doc comment.
var alpha_eff: f32;
if (is_biaxial && material.has_beta_ray != 0u) {
alpha_eff = pleochroic_channel_alpha_biaxial(
alpha_o_hoisted[k], alpha_beta_hoisted[k], alpha_e_hoisted[k], c_axis, current_plane_normal, wave_dir_at_bounce, eigen_a, eigen_b, stokes[k],
);
} else {
alpha_eff = pleochroic_channel_alpha(
alpha_o_hoisted[k], alpha_e_hoisted[k], c_axis, current_plane_normal, wave_dir_at_bounce, eigen_a, eigen_b, stokes[k],
);
}
// P1 (absorption path scale): mirrors
// optics::raytracer::absorption::apply_absorption's own
// `path_len * ctx.material.absorption_path_scale` multiply exactly.
let scaled_hit_t = hit.t * material.absorption_path_scale;
let trans_factor = exp(-alpha_eff * scaled_hit_t);
stokes[k] = stokes[k] * trans_factor;
}
}
}
// P2: angle of incidence measured against the WAVE NORMAL `k`
// (`wave_dir_at_bounce`), not the Poynting/energy direction `S` -- see this
// block's own design note above.
let cos_i = clamp(dot(-wave_dir_at_bounce, normal), 0.0, 1.0);
let sin_i = sqrt(max(fma(-cos_i, cos_i, 1.0), 0.0));
// Phase 3: theta_c fixed-point iteration (optics::raytracer::theta_c_for_bounce)
// plus the per-channel ordinary/effective-extraordinary index pair
// (optics::raytracer::per_channel_uniaxial_indices, called once per channel via
// per_channel_uniaxial_index -- see transport_physics.wgsl). For a cubic
// material (is_anisotropic == false) this reduces to n_eff_ch[k] == n_o_ch[k]
// for every k, exactly Phase 2's isotropic-only computation. For a biaxial
// material this is a DEAD (unused) computation -- see transport_physics.wgsl's
// `theta_c_for_bounce` doc comment. P2: fed `wave_dir_at_bounce` (`k`), not
// `current_dir` (`S`).
let theta_c = theta_c_for_bounce(normal, wave_dir_at_bounce, cos_i, inside_gem, is_anisotropic, c_axis, n_o_hero_seed, birefringence_delta);
// R3-GPU: `n_o_ch[k]` is exactly `n_o_hoisted[k]` (both `dispersion_evaluate` on
// the same `lambdas[k]`, hoisted above) -- only the theta_c-dependent
// `n_eff_k` half of `per_channel_uniaxial_index` still needs to run every
// bounce, so that's all this loop does now, reproducing that function's own
// `n_e_k`/`effective_extraordinary_index` body (see transport_physics.wgsl) on
// the hoisted `n_o_hoisted[k]` rather than calling the full function (which
// would redundantly re-evaluate dispersion) -- bit-identical either way.
//
// P3 (extraordinary-ray dispersion GPU port): `n_e_k` mirrors
// `optics::materials::GemMaterial::extraordinary_index_at` exactly -- a genuine
// wavelength-dependent evaluation of the material's OWN extraordinary-ray curve
// when `material.has_extraordinary_dispersion != 0` (Quartz/Amethyst/Citrine),
// falling back to the constant-offset `n_o_hoisted[k] + birefringence_delta`
// approximation every other uniaxial built-in still uses -- see
// `renderer::buffers::GpuGemMaterial`'s own doc comment.
var n_eff_ch: array<f32, 8>;
// P2 full uniaxial Fresnel (Lekner 1991): the RAW (not theta_c-projected)
// per-channel extraordinary index -- exactly `optics::materials::GemMaterial::
// extraordinary_index_at(lambdas[k], n_o_hoisted[k])` -- hoisted here (this
// loop already computes it as `n_e_k` on the way to `n_eff_k`) so the uniaxial
// entry/internal dispatch arms below can read it directly instead of
// recomputing the same dispersion evaluation a second time.
var n_e_raw_ch: array<f32, 8>;
for (var k: u32 = 0u; k < NUM_CHANNELS; k = k + 1u) {
var n_e_k: f32;
if (material.has_extraordinary_dispersion != 0u) {
n_e_k = extraordinary_dispersion_evaluate(material.extraordinary_model_type, material.extraordinary_param_a, material.extraordinary_param_b, lambdas[k]);
} else {
n_e_k = n_o_hoisted[k] + birefringence_delta;
}
n_e_raw_ch[k] = n_e_k;
var n_eff_k = n_o_hoisted[k];
if (is_anisotropic) {
n_eff_k = effective_extraordinary_index(n_o_hoisted[k], n_e_k, theta_c);
}
n_eff_ch[k] = n_eff_k;
}
let n_o_hero = n_o_hoisted[0];
let n_e_hero = n_o_hero + birefringence_delta;
// P2 full uniaxial Fresnel (Lekner 1991): the shared local frame + optic-axis
// direction cosines `uniaxial_fresnel`'s closed-form solver needs, built ONCE
// per bounce -- mirrors `BounceRefractionGeometry::uniaxial_frame`'s own "Some
// only when genuinely uniaxial" contract; harmless to build unconditionally
// here (cheap, pure function of already-hoisted inputs) since every use below
// is itself gated on `is_anisotropic && !is_biaxial`. `entry_incidence_frame`
// (performance, requirement 7: shared across all 8 channels of THIS bounce,
// `n1 == 1.0` fixed for an air->crystal entry) is likewise built once here even
// though only the entry arm below consumes it -- see that type's own doc
// comment.
let uframe = uniaxial_frame_build(wave_dir_at_bounce, normal, c_axis, cos_i, sin_i);
let uinc_frame = entry_incidence_frame(1.0, uframe);
// Mirrors `BounceRefractionGeometry::uniaxial_frame`'s own `Some` condition
// (`is_anisotropic && !is_biaxial`) exactly -- `apply_tir_bounce`'s CPU
// uniaxial branch gates on exactly this, with NO further degenerate-axis
// check (see below), so this flag mirrors it bit-for-bit including that
// narrow gap.
let uniaxial_active = is_anisotropic && !is_biaxial;
// Degenerate wave-normal-parallel-to-optic-axis limit -- mirrors
// `apply_partial_fresnel_bounce`'s identical `k_hat.cross(c_axis).
// length_squared()` guard exactly (see that function's own doc comment for why
// falling through to the existing scalar-at-`n_o` machinery is exact, not an
// approximation, at this exact limit). Used ONLY by the entry/general-internal
// dispatch arms below -- `apply_tir_bounce`'s own CPU uniaxial branch has NO
// equivalent guard (a pre-existing narrow gap this port deliberately mirrors
// rather than silently patches on one side only; see this task's own report).
let uniaxial_nondegenerate = uniaxial_active
&& dot(cross(wave_dir_at_bounce, c_axis), cross(wave_dir_at_bounce, c_axis)) > 1e-6;
// Phase 4: optics::raytracer::{hero_biaxial_wave_dirs, per_channel_biaxial_indices}.
// "mode A" is the faster (lower-index) root, "mode B" the slower -- see this
// file's header comment. Zero-initialized and never consulted downstream unless
// `is_biaxial` guards the read, exactly mirroring the CPU arrays' own
// "computed unconditionally, only ever POPULATED when is_biaxial" contract.
var wave_dir_a_hero = vec3<f32>(0.0, 0.0, 0.0);
var wave_dir_b_hero = vec3<f32>(0.0, 0.0, 0.0);
var n_biax_a_ch: array<f32, 8>;
var n_biax_b_ch: array<f32, 8>;
var n_alpha_ch: array<f32, 8>;
var n_gamma_ch: array<f32, 8>;
for (var k: u32 = 0u; k < NUM_CHANNELS; k = k + 1u) {
n_biax_a_ch[k] = 0.0;
n_biax_b_ch[k] = 0.0;
n_alpha_ch[k] = 0.0;
n_gamma_ch[k] = 0.0;
}
if (is_biaxial) {
if (inside_gem) {
// P2: both biaxial modes walk off (see optics::raytracer::refraction's
// hero_biaxial_wave_dirs doc comment) -- k, not S.
wave_dir_a_hero = wave_dir_at_bounce;
wave_dir_b_hero = wave_dir_at_bounce;
} else {
let res_a = biaxial_resolve_entry_mode(n_alpha_hero, n_beta_hero, n_gamma_hero, biax_ax0, biax_ax1, biax_ax2, wave_dir_at_bounce, normal, cos_i, n_o_hero_seed, false);
let res_b = biaxial_resolve_entry_mode(n_alpha_hero, n_beta_hero, n_gamma_hero, biax_ax0, biax_ax1, biax_ax2, wave_dir_at_bounce, normal, cos_i, n_o_hero_seed, true);
wave_dir_a_hero = res_a.wave_dir;
wave_dir_b_hero = res_b.wave_dir;
}
for (var k: u32 = 0u; k < NUM_CHANNELS; k = k + 1u) {
// optics::materials::GemMaterial::biaxial_indicatrix's per-channel
// n_beta/n_alpha/n_gamma -- n_beta_k is bit-identical to n_o_hoisted[k]
// (both `dispersion.evaluate(lambdas[k])`, same inputs), so reused
// rather than recomputed.
let n_beta_k = n_o_hoisted[k];
let n_alpha_k = n_beta_k - material.dispersion.biaxial_delta_beta_alpha;
let n_gamma_k = n_alpha_k + birefringence_delta;
n_alpha_ch[k] = n_alpha_k;
n_gamma_ch[k] = n_gamma_k;
let ni_a = biaxial_wave_indices(n_alpha_k, n_beta_k, n_gamma_k, biax_ax0, biax_ax1, biax_ax2, wave_dir_a_hero);
let ni_b = biaxial_wave_indices(n_alpha_k, n_beta_k, n_gamma_k, biax_ax0, biax_ax1, biax_ax2, wave_dir_b_hero);
n_biax_a_ch[k] = ni_a.y;
n_biax_b_ch[k] = ni_b.x;
}
}
let n_biax_a_hero = n_biax_a_ch[0];
let n_biax_b_hero = n_biax_b_ch[0];
// Fix 3c (CPU side) / Phase 3/4: while inside an anisotropic crystal, the medium
// index this ray is currently in is mode A or mode B depending on which
// eigenmode `is_extraordinary` selected at the most recent entry; outside the
// crystal (or for an isotropic material) it is always the mode-B array, which
// for a cubic material equals n_o_hoisted exactly (see per_channel_uniaxial_index).
// `is_biaxial` selects which pair of arrays ("mode A"/"mode B") is consulted --
// the biaxial ones computed just above, or the uniaxial n_o_hoisted/n_eff_ch pair.
let use_mode_a_medium = is_anisotropic && inside_gem && !is_extraordinary;
var n_medium_ch: array<f32, 8>;
for (var k: u32 = 0u; k < NUM_CHANNELS; k = k + 1u) {
var mode_a_k: f32;
var mode_b_k: f32;
if (is_biaxial) {
mode_a_k = n_biax_a_ch[k];
mode_b_k = n_biax_b_ch[k];
} else {
mode_a_k = n_o_hoisted[k];
mode_b_k = n_eff_ch[k];
}
n_medium_ch[k] = select(mode_b_k, mode_a_k, use_mode_a_medium);
}
let n_medium_hero = n_medium_ch[0];
let n1 = select(1.0, n_medium_hero, inside_gem);
let n2 = select(n_medium_hero, 1.0, inside_gem);
let eta = n1 / n2;
let sin2_t = eta * eta * fma(-cos_i, cos_i, 1.0);
// Task 2 GPU port (frosted girdle finish): which specular/diffuse treatment
// this facet gets. Bounds-checked (mirrors optics::raytracer::trace_spectral_ray_inner's
// `facet_finishes.get(hit_rec.facet_idx).copied().unwrap_or_default()`): an
// index past the end of a shorter-than-`planes` buffer defaults to
// `facet_finish::POLISHED`, exactly as the CPU's `Default` does.
// `select`'s two value arguments are BOTH evaluated in WGSL (no short-circuit),
// so an `if` guard is used instead of `select` here -- indexing
// `facet_finishes[hit.facet_idx]` unconditionally would attempt an
// out-of-bounds storage-buffer read whenever `facet_idx >= arrayLength(...)`
// (harmless per WGSL's robustness guarantees, but its RESULT is
// implementation-defined, so this must never be allowed to influence `finish`).
var finish: u32 = 0u;
if (hit.facet_idx < arrayLength(&facet_finishes)) {
finish = facet_finishes[hit.facet_idx];
}
// Captured before any of the three bounce-dispatch arms below run -- see
// optics::raytracer::dispatch_bounce's doc comment for why `was_internal_reflection`
// is always computed against the PRE-bounce `inside_gem`. The TIR/reflect/refract
// arms below never need this explicitly (they either never touch `inside_gem` at
// all, or only flip it once at their very end, after which nothing else in that
// arm reads it this bounce) -- only the frosted arm, which can update `inside_gem`
// itself, needs the pre-bounce value spelled out separately.
let pre_bounce_inside_gem = inside_gem;
if (finish == FACET_FINISH_FROSTED) {
// optics::raytracer::apply_frosted_bounce (via dispatch_bounce): the
// REPLACEMENT for the TIR/partial-reflect/refract dispatch below, not an
// addition to it -- see transport_physics.wgsl's own doc comment for the
// achromatic-by-design physics this must preserve exactly.
let fb = apply_frosted_bounce(
is_anisotropic, sin2_t, n1, n2, cos_i, normal, inside_gem, is_extraordinary,
seed0, bounce, &stokes, &path_pdf,
);
current_origin = hit_point + fb.new_dir * RAY_EPS;
current_dir = fb.new_dir;
// P2: a frosted (diffuse) bounce already depolarizes -- k collapses to S,
// mirroring optics::raytracer::transport::dispatch_bounce's identical
// treatment of FacetFinish::Frosted.
current_k = fb.new_dir;
inside_gem = fb.new_inside_gem != 0u;
if (fb.has_extraordinary_update != 0u) {
is_extraordinary = fb.extraordinary_update != 0u;
}
// Task 1: same was_internal_reflection formula optics::raytracer::dispatch_bounce
// uses for every bounce kind -- true for the frosted TIR-forced and reflect
// arms (inside_gem unchanged, no extraordinary update reported), false for
// the transmit arm (inside_gem always flips there).
let was_internal_reflection = pre_bounce_inside_gem
&& (fb.has_extraordinary_update == 0u)
&& (inside_gem == pre_bounce_inside_gem);
if (is_anisotropic && was_internal_reflection) {
// No stokes/path_pdf scaling here -- see internal_mode_coupling_draw's
// doc comment: this is a RELABELING of which eigenmode governs the
// NEXT bounce, not a SPLIT into two rays, so the matching unbiased
// scale factor is 1.0 (no-op) -- the same 1.0 conclusion the entry
// split itself reaches, for a different reason (see this file's header
// comment / `apply_refract_bounce`'s doc comment on the CPU side).
is_extraordinary = internal_mode_coupling_draw(
c_axis, is_biaxial, current_plane_normal, fb.new_dir,
stokes[0].x, stokes[0].y, stokes[0].z, false, 0.0, seed0, bounce,
);
}
} else if (sin2_t > 1.0) {
// Hero forced TIR (probability 1, no pdf division needed).
// P2 full uniaxial Fresnel (Lekner 1991): the closed-form o<->e-coupled
// TIR reflectance -- mirrors
// `optics::raytracer::refraction::apply_tir_bounce`'s own uniaxial branch
// bit-for-bit (including that function's own lack of a degenerate-axis
// guard, see `uniaxial_active`'s doc comment above -- this dispatch stays
// an EXACT mirror of the CPU function it replaces, not an independently
// "improved" one). `tir_exact_p_o` is the hero channel's own
// `R_o/(R_o+R_e)`, fed to `internal_mode_coupling_draw` below exactly as
// `apply_tir_bounce`'s own `exact_p_o` return value is on the CPU side.
var tir_exact_p_o: f32 = 0.5;
if (uniaxial_active) {
for (var k: u32 = 0u; k < NUM_CHANNELS; k = k + 1u) {
let n_inc_k = select(1.0, n_medium_ch[k], inside_gem);
let sol = internal_solve(n_inc_k, n_o_hoisted[k], n_e_raw_ch[k], c_axis, uframe, !is_extraordinary);
let flux_inc = max(sol.flux_inc, 1e-12);
let ro_pow = cplx_norm_sqr(sol.r_o) * sol.flux_ro;
let re_pow = cplx_norm_sqr(sol.r_e) * sol.flux_re;
let r_total_k = min((ro_pow + re_pow) / flux_inc, 1.0);
stokes[k] = stokes[k] * r_total_k;
let n2k_dbg = select(n_medium_ch[k], 1.0, inside_gem);
let etak_dbg = n_inc_k / n2k_dbg;
let sin2_t_k = etak_dbg * etak_dbg * fma(-cos_i, cos_i, 1.0);
if (sin2_t_k <= 1.0) {
path_pdf[k] = path_pdf[k] * clamp(r_total_k, 1e-4, 1.0 - 1e-4);
}
if (k == 0u) {
tir_exact_p_o = ro_pow / max(ro_pow + re_pow, 1e-12);
}
}
} else {
for (var k: u32 = 0u; k < NUM_CHANNELS; k = k + 1u) {
let n1k = select(1.0, n_medium_ch[k], inside_gem);
let n2k = select(n_medium_ch[k], 1.0, inside_gem);
let etak = n1k / n2k;
let sin2_t_k = etak * etak * fma(-cos_i, cos_i, 1.0);
if (sin2_t_k > 1.0) {
let delta_k = tir_phase_delta(n1k, cos_i, sin_i);
stokes[k] = mueller_tir_retardation(delta_k) * stokes[k];
} else {
let cos_t_k = sqrt(max(1.0 - sin2_t_k, 0.0));
let r_s_k = fma(n2k, -cos_t_k, n1k * cos_i) / fma(n2k, cos_t_k, n1k * cos_i);
let r_p_k = fma(n1k, -cos_t_k, n2k * cos_i) / fma(n1k, cos_t_k, n2k * cos_i);
stokes[k] = mueller_fresnel_reflection(r_s_k, r_p_k) * stokes[k];
let r_unpol_k = clamp(0.5 * fma(r_p_k, r_p_k, r_s_k * r_s_k), R_UNPOL_MIN, R_UNPOL_MAX);
path_pdf[k] = path_pdf[k] * r_unpol_k;
}
}
}
// P2: reflects the WAVE NORMAL `k` (not `S`), then re-derives `S'` for the
// reflected `k'` via `poynting_dir_for_mode` -- see this file's own design
// note above `wave_dir_at_bounce`.
let k_prime = wave_dir_at_bounce - 2.0 * dot(wave_dir_at_bounce, normal) * normal;
let s_prime = poynting_dir_for_mode(
is_anisotropic, is_biaxial, inside_gem, is_extraordinary, k_prime, c_axis,
n_o_hero, n_e_hero, n_alpha_hero, n_beta_hero, n_gamma_hero, biax_ax0, biax_ax1, biax_ax2,
);
current_origin = hit_point + s_prime * RAY_EPS;
current_dir = s_prime;
current_k = k_prime;
// Task 1: TIR is always an internal reflection (`n1 > n2` for the hero
// channel implies `inside_gem`, exactly as on the CPU side).
if (is_anisotropic) {
// No stokes/path_pdf scaling -- relabeling, not a split; see
// internal_mode_coupling_draw's doc comment.
is_extraordinary = internal_mode_coupling_draw(
c_axis, is_biaxial, current_plane_normal, k_prime,
stokes[0].x, stokes[0].y, stokes[0].z, uniaxial_active, tir_exact_p_o, seed0, bounce,
);
}
} else if (uniaxial_nondegenerate) {
// P2 full uniaxial Fresnel (Lekner 1991): entry AND general (non-forced)
// internal/exit dispatch, in full, self-contained -- mirrors
// optics::raytracer::refraction::apply_partial_fresnel_bounce's own early
// returns to apply_uniaxial_entry_bounce/apply_uniaxial_internal_bounce.
// `n_e_hero`/`n_e_ch` used elsewhere in this kernel for WALK-OFF/DIRECTION
// purposes are the constant-offset `n_o + birefringence_delta`
// approximation (matching `BounceRefractionGeometry::n_e_hero` on the CPU
// side); the closed-form FRESNEL solve below instead needs the material's
// own genuine extraordinary-ray dispersion curve where one exists (`n_e_
// raw_ch`/`n_e_hero_solve`, matching `GemMaterial::extraordinary_index_at`
// -- see `apply_uniaxial_entry_bounce`'s and `apply_uniaxial_internal_
// bounce`'s own local `n_e_hero` shadowing on the CPU side for why these
// are deliberately two different values, not a duplicate).
let n_e_hero_solve = n_e_raw_ch[0];
if (!inside_gem) {
// ---- ENTRY: mirrors apply_uniaxial_entry_bounce ----
let pair_hero = entry_solve_pair_with_incidence(uinc_frame, 1.0, n_o_hero, n_e_hero_solve, c_axis, uframe);
let sol_s_hero = pair_hero.s_sol;
let sol_p_hero = pair_hero.p_sol;
let reflect_mueller_hero = jones_to_mueller(sol_s_hero.r_s, sol_p_hero.r_s, sol_s_hero.r_p, sol_p_hero.r_p);
let entry_r_branch = clamp(reflect_mueller_hero[0][0], R_UNPOL_SELECT_MIN, R_UNPOL_SELECT_MAX);
let entry_inc_flux = 0.5 * uframe.cos_i;
let p_o_raw = mode_power(sol_s_hero.t_o, sol_p_hero.t_o, sol_s_hero.flux_o, entry_inc_flux, stokes[0]);
let p_e_raw = mode_power(sol_s_hero.t_e, sol_p_hero.t_e, sol_s_hero.flux_e, entry_inc_flux, stokes[0]);
let p_o_hero = clamp(p_o_raw / max(p_o_raw + p_e_raw, 1e-12), R_UNPOL_SELECT_MIN, R_UNPOL_SELECT_MAX);
let entry_mode_split_rand = f32(hash_u32(seed0 ^ hash_u32(bounce ^ BIREFRINGENT_SPLIT_STREAM))) / 4294967295.0;
let use_extraordinary_entry = entry_mode_split_rand < (1.0 - p_o_hero);
var p_mode_hero_frac = p_o_hero;
if (use_extraordinary_entry) {
p_mode_hero_frac = 1.0 - p_o_hero;
}
let n2_hero_dir_e = select(n_o_hero, n2, use_extraordinary_entry);
let eta_dir_e = n1 / n2_hero_dir_e;
let sin2_t_dir_e = min(eta_dir_e * eta_dir_e * fma(-cos_i, cos_i, 1.0), 1.0);
let cos_t_dir_e = sqrt(max(1.0 - sin2_t_dir_e, 0.0));
let refr_wave_dir_e = normalize(eta_dir_e * wave_dir_at_bounce + fma(eta_dir_e, cos_i, -cos_t_dir_e) * normal);
var final_refr_dir_e = refr_wave_dir_e;
if (use_extraordinary_entry) {
final_refr_dir_e = extraordinary_poynting_dir(refr_wave_dir_e, c_axis, n_o_hero, n_e_hero);
}
let p_axis_transmitted = cross(refr_wave_dir_e, uframe.s_axis);
let entry_rng_bounce = f32(hash_u32(seed0 ^ hash_u32(bounce ^ FRESNEL_BRANCH_STREAM))) / 4294967295.0;
if (entry_rng_bounce < entry_r_branch) {
for (var k: u32 = 0u; k < NUM_CHANNELS; k = k + 1u) {
let pair_k = entry_solve_pair_with_incidence(uinc_frame, 1.0, n_o_hoisted[k], n_e_raw_ch[k], c_axis, uframe);
let mueller_k = jones_to_mueller(pair_k.s_sol.r_s, pair_k.p_sol.r_s, pair_k.s_sol.r_p, pair_k.p_sol.r_p);
let r_unpol_k = clamp(mueller_k[0][0], 1e-4, 1.0 - 1e-4);
stokes[k] = (mueller_k * stokes[k]) * (1.0 / entry_r_branch);
path_pdf[k] = path_pdf[k] * r_unpol_k;
}
let new_k_e = wave_dir_at_bounce - 2.0 * dot(wave_dir_at_bounce, normal) * normal;
current_origin = hit_point + new_k_e * RAY_EPS;
current_dir = new_k_e;
current_k = new_k_e;
} else {
// P6: this is an INTERIOR dispersive event (air->crystal entry) --
// mismatch keeps `path_pdf` accumulating (via this event's own
// `t_unpol_k`) and narrows `compat` below, no split -- see this
// file's own header comment and
// `apply_uniaxial_entry_transmit_channels`'s CPU-side doc comment.
var entry_dirs: array<vec3<f32>, 8>;
var entry_dirs_valid: array<bool, 8>;
var entry_hero_match: array<bool, 8>;
for (var k: u32 = 0u; k < NUM_CHANNELS; k = k + 1u) {
let n_o_k = n_o_hoisted[k];
var n2k_dir = n_o_k;
if (use_extraordinary_entry) {
n2k_dir = n_eff_ch[k];
}
let eta_dir_k = n1 / n2k_dir;
let sin2_t_k = eta_dir_k * eta_dir_k * fma(-cos_i, cos_i, 1.0);
if (sin2_t_k > 1.0) {
stokes[k] = stokes[k] * 0.0;
path_pdf[k] = 0.0;
continue;
}
let cos_t_k = sqrt(max(1.0 - sin2_t_k, 0.0));
let refr_wave_dir_k = normalize(eta_dir_k * wave_dir_at_bounce + fma(eta_dir_k, cos_i, -cos_t_k) * normal);
var final_dir_k = refr_wave_dir_k;
if (use_extraordinary_entry) {
let n_e_eff_k = n_o_k + birefringence_delta;
final_dir_k = extraordinary_poynting_dir(refr_wave_dir_k, c_axis, n_o_k, n_e_eff_k);
}
entry_dirs[k] = final_dir_k;
entry_dirs_valid[k] = true;
let direction_matches_entry = dot(final_dir_k, final_refr_dir_e) >= DIRECTION_MATCH_COS_TOL;
entry_hero_match[k] = direction_matches_entry;
if (!direction_matches_entry) {
stokes[k] = stokes[k] * 0.0;
// P6: technique k stays a live member of every compatible
// channel's MIS family -- keep its own entry-transmit
// density accumulating (the SAME t_unpol_k the matching
// branch below folds in); only its radiance ends here.
let pair_k_mm = entry_solve_pair_with_incidence(uinc_frame, 1.0, n_o_k, n_e_raw_ch[k], c_axis, uframe);
let mueller_k_mm = jones_to_mueller(pair_k_mm.s_sol.r_s, pair_k_mm.p_sol.r_s, pair_k_mm.s_sol.r_p, pair_k_mm.p_sol.r_p);
let t_unpol_k_mm = clamp(1.0 - mueller_k_mm[0][0], 1e-4, 1.0 - 1e-4);
path_pdf[k] = path_pdf[k] * t_unpol_k_mm;
continue;
}
let pair_k = entry_solve_pair_with_incidence(uinc_frame, 1.0, n_o_k, n_e_raw_ch[k], c_axis, uframe);
var t_s_row = pair_k.s_sol.t_o;
var t_p_row = pair_k.p_sol.t_o;
var flux_k = pair_k.s_sol.flux_o;
var mode_dir = pair_k.s_sol.o_hat;
if (use_extraordinary_entry) {
t_s_row = pair_k.s_sol.t_e;
t_p_row = pair_k.p_sol.t_e;
flux_k = pair_k.s_sol.flux_e;
mode_dir = pair_k.s_sol.e_hat;
}
let p_mode_k = mode_power(t_s_row, t_p_row, flux_k, entry_inc_flux, stokes[k]);
let deposit_i = p_mode_k / p_mode_hero_frac;
let azimuth = azimuth2_in_frame(mode_dir, uframe.s_axis, p_axis_transmitted);
stokes[k] = vec4<f32>(deposit_i, deposit_i * azimuth.x, deposit_i * azimuth.y, 0.0) * (1.0 / (1.0 - entry_r_branch));
let mueller_k = jones_to_mueller(pair_k.s_sol.r_s, pair_k.p_sol.r_s, pair_k.s_sol.r_p, pair_k.p_sol.r_p);
let t_unpol_k = clamp(1.0 - mueller_k[0][0], 1e-4, 1.0 - 1e-4);
path_pdf[k] = path_pdf[k] * t_unpol_k;
}
narrow_compat(&compat, entry_dirs, entry_dirs_valid, entry_hero_match);
current_origin = hit_point + final_refr_dir_e * RAY_EPS;
current_dir = final_refr_dir_e;
current_k = refr_wave_dir_e;
inside_gem = true;
is_extraordinary = use_extraordinary_entry;
}
} else {
// ---- INTERNAL: mirrors apply_uniaxial_internal_bounce ----
let sol_hero_i = internal_solve(n1, n_o_hero, n_e_hero_solve, c_axis, uframe, !is_extraordinary);
let flux_inc_hero_i = max(sol_hero_i.flux_inc, 1e-12);
let ro_pow_hero = cplx_norm_sqr(sol_hero_i.r_o) * sol_hero_i.flux_ro;
let re_pow_hero = cplx_norm_sqr(sol_hero_i.r_e) * sol_hero_i.flux_re;
let internal_r_branch = clamp((ro_pow_hero + re_pow_hero) / flux_inc_hero_i, R_UNPOL_SELECT_MIN, R_UNPOL_SELECT_MAX);
let internal_p_o_exact = ro_pow_hero / max(ro_pow_hero + re_pow_hero, 1e-12);
let internal_rng_bounce = f32(hash_u32(seed0 ^ hash_u32(bounce ^ FRESNEL_BRANCH_STREAM))) / 4294967295.0;
if (internal_rng_bounce < internal_r_branch) {
for (var k: u32 = 0u; k < NUM_CHANNELS; k = k + 1u) {
var sol = sol_hero_i;
if (k != 0u) {
sol = internal_solve(n_medium_ch[k], n_o_hoisted[k], n_e_raw_ch[k], c_axis, uframe, !is_extraordinary);
}
let flux_inc_k = max(sol.flux_inc, 1e-12);
let r_total_k = min((cplx_norm_sqr(sol.r_e) * sol.flux_re + cplx_norm_sqr(sol.r_o) * sol.flux_ro) / flux_inc_k, 1.0);
stokes[k] = stokes[k] * (r_total_k / internal_r_branch);
path_pdf[k] = path_pdf[k] * clamp(r_total_k, 1e-4, 1.0 - 1e-4);
}
let new_k_i = wave_dir_at_bounce - 2.0 * dot(wave_dir_at_bounce, normal) * normal;
current_origin = hit_point + new_k_i * RAY_EPS;
current_dir = new_k_i;
current_k = new_k_i;
is_extraordinary = internal_mode_coupling_draw(
c_axis, is_biaxial, current_plane_normal, new_k_i,
stokes[0].x, stokes[0].y, stokes[0].z, true, internal_p_o_exact, seed0, bounce,
);
} else {
let eta_dir_i = n1 / n2;
let sin2_t_dir_i = min(eta_dir_i * eta_dir_i * fma(-cos_i, cos_i, 1.0), 1.0);
let cos_t_dir_i = sqrt(max(1.0 - sin2_t_dir_i, 0.0));
let hero_refr_dir_i = normalize(eta_dir_i * wave_dir_at_bounce + fma(eta_dir_i, cos_i, -cos_t_dir_i) * normal);
// P6: this IS the uniaxial exact EXIT event (crystal -> air) --
// mismatch resolves its own split (point 1) and keeps `path_pdf`
// accumulating via this event's own exit factor (point 2), but
// never narrows `compat` -- see
// `apply_uniaxial_internal_transmit_channels`'s CPU-side doc
// comment and this file's own header comment.
for (var k: u32 = 0u; k < NUM_CHANNELS; k = k + 1u) {
// P6: captured before any mutation below, so the split branch
// (which needs the ORIGINAL incident value after `stokes[k]`
// has already been zeroed by the unchanged chromatic-
// termination bookkeeping) can still compute k's own
// transmission.
let original_stokes_i = stokes[k].x;
let n_inc_k = n_medium_ch[k];
let eta_dir_k = n_inc_k;
let sin2_t_k = eta_dir_k * eta_dir_k * fma(-cos_i, cos_i, 1.0);
if (sin2_t_k > 1.0) {
stokes[k] = stokes[k] * 0.0;
path_pdf[k] = 0.0;
continue;
}
let cos_t_k = sqrt(max(1.0 - sin2_t_k, 0.0));
let refr_wave_dir_k = normalize(eta_dir_k * wave_dir_at_bounce + fma(eta_dir_k, cos_i, -cos_t_k) * normal);
let direction_matches_i = dot(refr_wave_dir_k, hero_refr_dir_i) >= DIRECTION_MATCH_COS_TOL;
// Performance: reuse the caller's already-solved hero channel
// instead of re-solving the identical boundary system -- needed
// by both branches below (matching directly; split via
// `compute_uniaxial_exit_transmission`).
var sol = sol_hero_i;
if (k != 0u) {
sol = internal_solve(n_inc_k, n_o_hoisted[k], n_e_raw_ch[k], c_axis, uframe, !is_extraordinary);
}
if (!direction_matches_i) {
// Chromatic termination -- UNCHANGED from pre-P6.
let prefix_path_pdf_k = path_pdf[k];
stokes[k] = stokes[k] * 0.0;
path_pdf[k] = 0.0;
let et_mm = compute_uniaxial_exit_transmission(sol, internal_r_branch, original_stokes_i);
let t_unpol_k_mm = clamp(et_mm.i_unit, 1e-4, 1.0 - 1e-4);
path_pdf[k] = prefix_path_pdf_k * t_unpol_k_mm;
if (original_stokes_i > 0.0) {
try_split_exit_channel(
&split_radiance, hit_point, k, lambdas[k], refr_wave_dir_k,
et_mm.transmitted.x, studio_key_dir, studio_fill_dir, studio_sin_lp,
);
}
continue;
}
// direction_matches: exactly the pre-P6 matching-branch formula.
let et = compute_uniaxial_exit_transmission(sol, internal_r_branch, original_stokes_i);
stokes[k] = et.transmitted;
let t_unpol_k = clamp(et.i_unit, 1e-4, 1.0 - 1e-4);
path_pdf[k] = path_pdf[k] * t_unpol_k;
}
current_origin = hit_point + hero_refr_dir_i * RAY_EPS;
current_dir = hero_refr_dir_i;
current_k = hero_refr_dir_i;
inside_gem = false;
}
}
} else {
// P1/P2: at an air->crystal entry into an anisotropic material, which
// eigenmode (ordinary/mode-A vs extraordinary/mode-B) this path's single
// geometric transmission event represents is decided HERE -- before the
// reflect-vs-transmit draw below -- rather than inside the REFRACT arm's
// own transmit branch as it used to be. See
// optics::raytracer::refraction::apply_partial_fresnel_bounce's doc
// comment (CPU side) for the full two-part rationale: (P1) the
// polarization-weighted selection (`entry_eigenmode_selection`) must be
// shared by both branches below -- a beam already aligned with the
// ordinary axis should be MORE likely to reflect at the ordinary index
// too, not just more likely to transmit as ordinary conditional on
// transmitting -- and (P2) the REFLECT branch's own Fresnel coefficients
// must be evaluated at the SAME mode's index the transmit branch uses, so
// `R + T == 1` for whichever mode this draw actually selects. A biaxial
// material has no uniaxial "ordinary" eigenmode to weight against, so both
// the selection and the index correction below are gated on `!is_biaxial`.
let entering_anisotropic = !inside_gem && is_anisotropic;
var entry_valid = 0u;
var entry_cos_2psi_o = 0.0;
var entry_sin_2psi_o = 0.0;
var p_o = 0.5;
if (entering_anisotropic && !is_biaxial) {
let sel = entry_eigenmode_selection(
c_axis, current_plane_normal, wave_dir_at_bounce,
stokes[0].x, stokes[0].y, stokes[0].z,
);
entry_valid = sel.valid;
entry_cos_2psi_o = sel.cos_2psi_o;
entry_sin_2psi_o = sel.sin_2psi_o;
if (sel.valid != 0u) {
p_o = sel.p_o;
}
}
let mode_split_rand = f32(hash_u32(seed0 ^ hash_u32(bounce ^ BIREFRINGENT_SPLIT_STREAM))) / 4294967295.0;
var use_extraordinary = is_extraordinary;
if (entering_anisotropic) {
// Must match optics::raytracer::refraction::apply_partial_fresnel_bounce's
// `mode_split_rand < (1.0 - p_o)` EXACTLY (same operator, same sense) --
// not just the same marginal probability. `rand < 1-p_o` and `rand >= p_o`
// both integrate to the same P(extraordinary) = 1-p_o, but they pick
// DIFFERENT (at p_o==0.5 exactly complementary) halves of the same [0,1)
// draw, so for an identical (seed0, bounce) -- and hence an identical
// `mode_split_rand` bit pattern on both CPU and GPU -- the two forms
// disagree on which mode a given path actually gets. This previously read
// `mode_split_rand >= p_o`, silently breaking CPU/GPU parity at every
// anisotropic entry (most visibly for unpolarized light, p_o==0.5: the old
// form picked extraordinary for rand>=0.5, the CPU/original convention
// picks it for rand<0.5 -- exact opposites, not just a relabeling).
use_extraordinary = mode_split_rand < (1.0 - p_o);
}
// The chosen mode's own doubled polarization azimuth, for the per-channel
// eigenmode projection in the REFRACT arm below -- meaningless unless
// `entry_valid != 0u`. Extraordinary is perpendicular to ordinary
// (psi_e == psi_o + 90deg), so its doubled azimuth is the negation of the
// ordinary one.
var entry_cos_2psi_x = entry_cos_2psi_o;
var entry_sin_2psi_x = entry_sin_2psi_o;
if (use_extraordinary) {
entry_cos_2psi_x = -entry_cos_2psi_o;
entry_sin_2psi_x = -entry_sin_2psi_o;
}
// P2: use the SELECTED mode's own index for n2 at THIS interface -- `n2`
// (mode B/extraordinary) unchanged when extraordinary is selected, not
// entering an anisotropic material at all, or biaxial (`n_o_hero` is a
// uniaxial-only quantity, not mode A), `n_o_hero` when the ordinary mode
// is selected instead. Bit-identical to the pre-P2 `sqrt(1.0 - sin2_t)` in
// every case this reduces to (`n2_decision == n2`): both `n1 / n2_decision`
// and `sin2_t` itself are pure functions of `n1`/`n2`/`cos_i`, so recomputing
// the ratio here rather than reusing the hero-level `sin2_t` produces the
// identical bits when the index is unchanged.
let n2_decision = select(n2, n_o_hero, entering_anisotropic && !is_biaxial && !use_extraordinary);
let eta_decision = n1 / n2_decision;
let cos_t = sqrt(max(1.0 - eta_decision * eta_decision * fma(-cos_i, cos_i, 1.0), 0.0));
let r_s = fma(n2_decision, -cos_t, n1 * cos_i) / fma(n2_decision, cos_t, n1 * cos_i);
let r_p = fma(n1, -cos_t, n2_decision * cos_i) / fma(n1, cos_t, n2_decision * cos_i);
let r_unpol_raw = 0.5 * fma(r_p, r_p, r_s * r_s);
// P5: [0.02, 0.98] here (distinct from the per-channel R_UNPOL_MIN/MAX
// used only to scale path_pdf, never to divide stokes) caps the
// `1/r_unpol`/`1/(1-r_unpol)` divisions below at 50x instead of 10,000x at
// grazing incidence -- still unbiased, just far less firefly-prone.
let r_unpol = clamp(r_unpol_raw, R_UNPOL_SELECT_MIN, R_UNPOL_SELECT_MAX);
let rng_bounce = f32(hash_u32(seed0 ^ hash_u32(bounce ^ FRESNEL_BRANCH_STREAM))) / 4294967295.0;
if (rng_bounce < r_unpol) {
// REFLECT: each channel applies its OWN Fresnel/TIR matrix, divided by
// the SAME hero selection probability r_unpol (the PDF-division
// coupling that makes GPU/CPU float divergence harmless -- see this
// file's header comment).
for (var k: u32 = 0u; k < NUM_CHANNELS; k = k + 1u) {
let n1k = select(1.0, n_medium_ch[k], inside_gem);
let n2k = select(n_medium_ch[k], 1.0, inside_gem);
let etak = n1k / n2k;
let sin2_t_k = etak * etak * fma(-cos_i, cos_i, 1.0);
var refl: mat4x4<f32>;
if (sin2_t_k > 1.0) {
let delta_k = tir_phase_delta(n1k, cos_i, sin_i);
refl = mueller_tir_retardation(delta_k);
} else {
let cos_t_k = sqrt(max(1.0 - sin2_t_k, 0.0));
let r_s_k = fma(n2k, -cos_t_k, n1k * cos_i) / fma(n2k, cos_t_k, n1k * cos_i);
let r_p_k = fma(n1k, -cos_t_k, n2k * cos_i) / fma(n1k, cos_t_k, n2k * cos_i);
let r_unpol_k = clamp(0.5 * fma(r_p_k, r_p_k, r_s_k * r_s_k), R_UNPOL_MIN, R_UNPOL_MAX);
path_pdf[k] = path_pdf[k] * r_unpol_k;
refl = mueller_fresnel_reflection(r_s_k, r_p_k);
}
stokes[k] = (refl * stokes[k]) * (1.0 / r_unpol);
}
// P2: reflects the WAVE NORMAL `k` (not `S`), then re-derives `S'` for
// the reflected `k'` via `poynting_dir_for_mode` -- see this file's own
// design note above `wave_dir_at_bounce`.
let k_prime = wave_dir_at_bounce - 2.0 * dot(wave_dir_at_bounce, normal) * normal;
let s_prime = poynting_dir_for_mode(
is_anisotropic, is_biaxial, inside_gem, is_extraordinary, k_prime, c_axis,
n_o_hero, n_e_hero, n_alpha_hero, n_beta_hero, n_gamma_hero, biax_ax0, biax_ax1, biax_ax2,
);
current_origin = hit_point + s_prime * RAY_EPS;
current_dir = s_prime;
current_k = k_prime;
// Task 1: this arm never changes `inside_gem` (unlike the refract arm
// below, which always flips it), so `inside_gem` here still holds its
// pre-bounce value -- an internal reflection iff it was already true.
if (is_anisotropic && inside_gem) {
// No stokes/path_pdf scaling -- relabeling, not a split; see
// internal_mode_coupling_draw's doc comment. This arm is only
// reached for a biaxial material or a degenerate-axis uniaxial one
// (the non-degenerate uniaxial case takes the closed-form
// `uniaxial_nondegenerate` branch above instead), so no exact p_o
// is available here -- falls back to the polarization-projection
// heuristic exactly as before.
is_extraordinary = internal_mode_coupling_draw(
c_axis, is_biaxial, current_plane_normal, k_prime,
stokes[0].x, stokes[0].y, stokes[0].z, false, 0.0, seed0, bounce,
);
}
} else {
// REFRACT -- optics::raytracer::apply_refract_bounce. P1/P2:
// `entering_anisotropic` and `use_extraordinary` are ALREADY resolved
// above (shared with the REFLECT arm's own r_unpol computation) -- no
// fresh draw here. On an air->crystal entry into an anisotropic
// material, the selected mode's own Fresnel transmittance is applied
// to that mode's PROJECTED Stokes state (see the per-channel loop
// below) at its own unscaled intensity: the selection probability is
// set to match this mode's own true physical energy fraction exactly,
// so multiplying by that fraction and dividing by the identical
// selection probability cancel -- not the naive `1/0.5` a disjoint
// split would need. For a cubic material (or any bounce that is not an
// anisotropic entry) `entering_anisotropic` is false and
// `use_extraordinary` keeps whatever `is_extraordinary` already was.
// Direction: the mode-A eigenmode uses n_mode_a and (uniaxial only) is
// never walked off; the mode-B eigenmode's ENERGY (Poynting) direction
// is displaced by the walk-off angle -- computed BEFORE the per-channel
// loop below so each companion channel's own hypothetical direction can
// be compared against this SAME hero-driven direction (Fix G / Part 2 on
// the CPU side). Phase 4: for a biaxial material entering the crystal,
// NEITHER mode is a plain constant-index Snell refraction -- BOTH modes
// walk off via `biaxial_mode_poynting_dir`, using `n_biax_a_hero`/
// `n_biax_b_hero` (the SAME looked-up scalars the per-channel loop's own
// `k == hero_idx` iteration uses) for self-consistency.
// P2: `refr_wave_dir` is the SNELL-REFRACTED WAVE NORMAL `k'` (fed from
// `wave_dir_at_bounce`, not `current_dir`/`S` -- see this file's
// own design note above `wave_dir_at_bounce`, and rule 4 in particular:
// "at exit into air, refract k (not S)"). Captured into `new_k` alongside
// the Poynting-converted `final_refr_dir` (`S'`), since the caller needs
// BOTH from here on.
var new_k: vec3<f32>;
var final_refr_dir: vec3<f32>;
if (entering_anisotropic && is_biaxial) {
let n2_hero_dir = select(n_biax_a_hero, n_biax_b_hero, use_extraordinary);
let eta_dir = n1 / n2_hero_dir;
let sin2_t_dir = min(eta_dir * eta_dir * fma(-cos_i, cos_i, 1.0), 1.0);
let cos_t_dir = sqrt(max(1.0 - sin2_t_dir, 0.0));
let refr_wave_dir = normalize(
eta_dir * wave_dir_at_bounce + fma(eta_dir, cos_i, -cos_t_dir) * normal,
);
new_k = refr_wave_dir;
final_refr_dir = biaxial_mode_poynting_dir(n_alpha_hero, n_beta_hero, n_gamma_hero, biax_ax0, biax_ax1, biax_ax2, refr_wave_dir, use_extraordinary);
} else {
let n2_hero_dir = select(n2, n_o_hero, entering_anisotropic && !use_extraordinary);
let eta_dir = n1 / n2_hero_dir;
let sin2_t_dir = min(eta_dir * eta_dir * fma(-cos_i, cos_i, 1.0), 1.0);
let cos_t_dir = sqrt(max(1.0 - sin2_t_dir, 0.0));
let refr_wave_dir = normalize(
eta_dir * wave_dir_at_bounce + fma(eta_dir, cos_i, -cos_t_dir) * normal,
);
new_k = refr_wave_dir;
final_refr_dir = refr_wave_dir;
if (entering_anisotropic && use_extraordinary) {
final_refr_dir = extraordinary_poynting_dir(refr_wave_dir, c_axis, n_o_hero, n_e_hero);
}
}
// P6: `is_exit_event` mirrors apply_refract_channel's own
// `inside_gem && !entering_anisotropic` -- since `entering_anisotropic`
// is only ever true while `!inside_gem`, this reduces to plain
// `inside_gem` (still the PRE-flip value here, unchanged until after
// this loop). Interior (entry) mismatches instead narrow `compat`,
// never split -- see this file's own header comment.
let is_exit_event = inside_gem;
var scalar_dirs: array<vec3<f32>, 8>;
var scalar_dirs_valid: array<bool, 8>;
var scalar_hero_match: array<bool, 8>;
for (var k: u32 = 0u; k < NUM_CHANNELS; k = k + 1u) {
// P6: captured before any mutation below, so the mismatch/split
// branch (which needs the ORIGINAL incident value after `stokes[k]`
// has already been zeroed by the unchanged chromatic-termination
// bookkeeping) can still compute k's own transmission.
let original_stokes_k = stokes[k];
let n1k = select(1.0, n_medium_ch[k], inside_gem);
var n2k = select(n_medium_ch[k], 1.0, inside_gem);
if (entering_anisotropic && !use_extraordinary) {
if (is_biaxial) {
n2k = n_biax_a_ch[k];
} else {
n2k = n_o_hoisted[k];
}
}
let ratio_k = n1k / n2k;
let sin2_t_k = (ratio_k * ratio_k) * fma(-cos_i, cos_i, 1.0);
if (sin2_t_k > 1.0) {
// Chromatic termination: channel k cannot transmit at this
// angle even though the hero-driven path did. Both the Stokes
// contribution AND the path_pdf are dropped to exactly 0 --
// never just down-weighted (see this file's header comment).
stokes[k] = stokes[k] * 0.0;
path_pdf[k] = 0.0;
continue;
}
let cos_t_k = sqrt(max(1.0 - sin2_t_k, 0.0));
let refr_wave_dir_k = normalize(
ratio_k * wave_dir_at_bounce + fma(ratio_k, cos_i, -cos_t_k) * normal,
);
// Fix G (Part 2) / the direction-match identity trap: channel k's own
// walk-off, using k's own per-channel indices, compared against the
// STORED `final_refr_dir` above -- never a second recomputation of
// the hero's own direction (which would be a few ULP different and
// chromatically self-terminate the hero channel against itself).
// Phase 4: channel k's own biaxial walk-off, using k's own
// (n_alpha_ch[k], n_beta_ch := n_o_hoisted[k], n_gamma_ch[k]) evaluated
// at k's own single-shot refracted wave direction -- the direct
// per-channel generalization of the uniaxial extraordinary_poynting_dir
// call below.
var final_dir_k = refr_wave_dir_k;
if (entering_anisotropic && is_biaxial) {
final_dir_k = biaxial_mode_poynting_dir(n_alpha_ch[k], n_o_hoisted[k], n_gamma_ch[k], biax_ax0, biax_ax1, biax_ax2, refr_wave_dir_k, use_extraordinary);
} else if (entering_anisotropic && use_extraordinary) {
let n_e_k = n_o_hoisted[k] + birefringence_delta;
final_dir_k = extraordinary_poynting_dir(refr_wave_dir_k, c_axis, n_o_hoisted[k], n_e_k);
}
scalar_dirs[k] = final_dir_k;
scalar_dirs_valid[k] = true;
let direction_matches = dot(final_dir_k, final_refr_dir) >= DIRECTION_MATCH_COS_TOL;
scalar_hero_match[k] = direction_matches;
if (direction_matches) {
let t_s_k = (2.0 * n1k * cos_i) / fma(n2k, cos_t_k, n1k * cos_i);
let t_p_k = (2.0 * n1k * cos_i) / fma(n1k, cos_t_k, n2k * cos_i);
let trans = mueller_fresnel_transmission(n1k, n2k, cos_i, cos_t_k, t_s_k, t_p_k);
// P1: project this channel's incident Stokes state onto the
// SELECTED eigenmode before transmission -- see
// entry_eigenmode_selection's doc comment for the full
// derivation. Fully linear along the mode's own axis (`w`/`V`
// zeroed), at this channel's OWN unscaled intensity (the
// mode-selection draw's probability matches this mode's true
// physical energy fraction exactly, so multiplying by that
// fraction and dividing by the identical selection probability
// cancel). `entry_valid == 0u` (biaxial, or negligible linear
// polarization) leaves `stokes[k]` unprojected, reducing
// exactly to the pre-P1 behaviour.
var incident_k = stokes[k];
if (entering_anisotropic && entry_valid != 0u) {
let i_k = stokes[k].x;
incident_k = vec4<f32>(i_k, i_k * entry_cos_2psi_x, i_k * entry_sin_2psi_x, 0.0);
}
// No `/ split_pdf` -- see the entering_anisotropic comment above.
stokes[k] = (trans * incident_k) * (1.0 / (1.0 - r_unpol));
let r_s_k = fma(n2k, -cos_t_k, n1k * cos_i) / fma(n2k, cos_t_k, n1k * cos_i);
let r_p_k = fma(n1k, -cos_t_k, n2k * cos_i) / fma(n1k, cos_t_k, n2k * cos_i);
let r_unpol_k = clamp(0.5 * fma(r_p_k, r_p_k, r_s_k * r_s_k), R_UNPOL_MIN, R_UNPOL_MAX);
// No `* split_pdf` -- scale-invariant under a uniform per-channel
// factor, was a pure no-op on the MIS weight; see refraction.rs.
path_pdf[k] = path_pdf[k] * (1.0 - r_unpol_k);
} else {
// P6: chromatic termination -- UNCHANGED radiance handling
// (mismatched channel loses its own Stokes contribution
// either way), but `path_pdf[k]` keeps accumulating this
// event's own transmit factor and, at an EXIT event, the
// channel additionally resolves its own transmitted radiance
// along its own refracted direction via
// `try_split_exit_channel` -- see this file's own header
// comment and `apply_refract_channel`'s CPU-side `else` branch.
let prefix_path_pdf_k = path_pdf[k];
stokes[k] = stokes[k] * 0.0;
path_pdf[k] = 0.0;
let azimuth_valid_here = entering_anisotropic && entry_valid != 0u;
let ct = compute_channel_transmission(
n1k, n2k, cos_i, cos_t_k, r_unpol, entering_anisotropic,
azimuth_valid_here, entry_cos_2psi_x, entry_sin_2psi_x, original_stokes_k,
);
path_pdf[k] = prefix_path_pdf_k * (1.0 - ct.r_unpol_k);
if (is_exit_event && original_stokes_k.x > 0.0) {
try_split_exit_channel(
&split_radiance, hit_point, k, lambdas[k], refr_wave_dir_k,
ct.transmitted.x, studio_key_dir, studio_fill_dir, studio_sin_lp,
);
}
}
}
// P6: an INTERIOR dispersive event (an entry into the gem) narrows
// every channel's MIS family; the exit event never does -- see
// `narrow_compat`'s own comment in `transport_physics.wgsl`.
if (!inside_gem) {
narrow_compat(&compat, scalar_dirs, scalar_dirs_valid, scalar_hero_match);
}
current_origin = hit_point + final_refr_dir * RAY_EPS;
current_dir = final_refr_dir;
current_k = new_k;
inside_gem = !inside_gem;
if (entering_anisotropic) {
is_extraordinary = use_extraordinary;
}
}
}
if (bounce > 4u) {
var max_intensity: f32 = 0.0;
for (var k: u32 = 0u; k < NUM_CHANNELS; k = k + 1u) {
max_intensity = max(max_intensity, max(stokes[k].x, 0.0));
}
let q = clamp(max_intensity, RR_FLOOR, 1.0);
let rr_rand = f32(hash_u32(seed0 ^ hash_u32(bounce ^ RUSSIAN_ROULETTE_STREAM))) / 4294967295.0;
if (rr_rand > q) {
break;
}
// P6 (Bias A fix): `split_radiance` rides along on the same `1/q` survival
// rescale as `stokes` -- see `apply_russian_roulette`'s CPU-side doc comment.
for (var k: u32 = 0u; k < NUM_CHANNELS; k = k + 1u) {
stokes[k] = stokes[k] * (1.0 / q);
split_radiance[k] = split_radiance[k] / q;
}
}
}
// P6: commit the staged exit-split contributions into the REAL `radiance` array
// ONLY if the shared/hero path itself reached its own environment lookup -- see
// `split_radiance`'s own comment above and `ExitSplitCtx::split_radiance`'s
// CPU-side doc comment for why this all-or-nothing gate is required. A no-op
// whenever nothing ever split, or every split channel's own probe declined.
if (path_escaped) {
for (var k: u32 = 0u; k < NUM_CHANNELS; k = k + 1u) {
radiance[k] = radiance[k] + split_radiance[k];
}
}
// P6: optics::raytracer::color::integrate_channels_to_xyz_families -- the
// per-channel-family generalisation of the pre-P6 single shared `spectral_mis_weight`
// above. Every channel that is alive at the exit (matching or split) contributes its
// own radiance, so each channel's balance-heuristic weight is normalised over
// exactly the techniques (hero choices) under which THAT channel would have stayed
// alive on this same geometric path -- its own family, `compat[k]` -- rather than
// the hero's. See `integrate_channels_to_xyz_families`'s own doc comment for why a
// single shared weight is biased once companions contribute, and why this reduces
// to it exactly whenever every family is the full set.
var xyz = vec3<f32>(0.0, 0.0, 0.0);
for (var k: u32 = 0u; k < NUM_CHANNELS; k = k + 1u) {
var family_pdf: f32 = 0.0;
for (var j: u32 = 0u; j < NUM_CHANNELS; j = j + 1u) {
if ((compat[k] & (1u << j)) != 0u) {
family_pdf = family_pdf + path_pdf[j];
}
}
// Same "should not happen" fallback as the pre-P6 shared `spectral_mis_weight`.
var weight_k: f32 = 1.0;
if (family_pdf > 1e-12) {
weight_k = f32(NUM_CHANNELS) * path_pdf[0] / family_pdf;
}
let cmf = cie_1931_cmf(lambdas[k]);
let weighted = radiance[k] * weight_k;
xyz = xyz + cmf * (weighted * NORM_FACTOR);
}
if (params.env_mode == 1u) {
// Fix 3: Bradford-LMS-space von Kries adaptation, not a raw XYZ scale -- see
// `apply_von_kries_white_balance`'s doc comment above and
// `optics::raytracer::apply_von_kries_white_balance` on the CPU side.
xyz = apply_von_kries_white_balance(xyz, params.white_balance);
}
// Mirrors `transport::trace_spectral_ray_inner`'s unconditional `.max(Vec3::ZERO)`
// on its white-balanced result (P6): the Bradford LMS round trip has negative
// off-diagonal entries and can push one component of a spectrally narrow sample
// slightly below zero (measured: Z = -4.9e-4 next to X = 1.75 on a diamond sample,
// where the CPU already reported exactly 0). Without this clamp the two engines
// disagree on such samples by that whole amount -- caught by the worker's
// `hybrid_gpu_and_cpu_split_sums_to_the_same_result_as_tracing_the_range_directly`
// as a 0.3 % Z-only drift of a 32-sample pixel sum. Applied outside the
// `env_mode` branch exactly as the CPU does; pre-white-balance `xyz` is a sum of
// non-negative terms, so the clamp is a no-op there.
xyz = max(xyz, vec3<f32>(0.0));
out_xyz[idx * 3u + 0u] = xyz.x;
out_xyz[idx * 3u + 1u] = xyz.y;
out_xyz[idx * 3u + 2u] = xyz.z;
// R4: `out_xyz` above is the only buffer a production dispatch
// (`renderer::gpu::frame::GpuFrameRenderer::accumulate`) ever reads back -- the four
// per-channel debug buffers below exist for Tier 2/spectral-debug self-tests only.
// Guarding their writes on `params.write_debug_buffers` (nonzero for every self-test,
// via `GpuTransportParams::new`'s default) lets a production dispatch bind tiny
// fixed-size dummy buffers for them instead of buffers sized like `out_xyz` -- 9x
// less write traffic and 9x more samples per chunk-budget dispatch. Every self-test's
// own behaviour is UNCHANGED: `write_debug_buffers != 0u` takes the exact same branch
// that ran unconditionally before this guard was added.
//
// P6: `out_compat` (one `compat[k]` mask per channel, same indexing as the other
// three) lets `estimator_check::run_spectral_debug`'s CPU-side self-consistency
// recombination call the real `integrate_channels_to_xyz_families` with the SAME
// per-channel family this dispatch itself used, instead of only the pre-P6
// `integrate_channels_to_xyz` (which would no longer agree with `out_xyz` whenever a
// trace ever narrowed a family).
if (params.write_debug_buffers != 0u) {
for (var k: u32 = 0u; k < NUM_CHANNELS; k = k + 1u) {
out_radiance[idx * 8u + k] = radiance[k];
out_lambdas[idx * 8u + k] = lambdas[k];
out_path_pdf[idx * 8u + k] = path_pdf[k];
out_compat[idx * 8u + k] = compat[k];
}
}
}