From 1ed95982ad2b336ebe40a538c47986275cdf58f2 Mon Sep 17 00:00:00 2001 From: Connor Johnstone Date: Fri, 16 Oct 2026 15:05:04 -0400 Subject: [PATCH] Added noise --- crates/resolved-cli/src/main.rs | 13 +- crates/resolved-cli/src/scene_file.rs | 19 +- crates/resolved/src/camera.rs | 63 ++++ crates/resolved/src/fidelity.rs | 5 +- crates/resolved/src/gpu/sensor.rs | 29 +- crates/resolved/src/render.rs | 29 +- crates/resolved/src/shaders/sensor.wgsl | 67 ++++- crates/resolved/tests/gpu/distortion.rs | 2 +- crates/resolved/tests/gpu/main.rs | 1 + crates/resolved/tests/gpu/noise.rs | 364 ++++++++++++++++++++++++ crates/resolved/tests/gpu/sensor.rs | 4 +- crates/resolved/tests/gpu/stars.rs | 3 + docs/optics.md | 101 ++++++- readme.md | 9 +- 14 files changed, 687 insertions(+), 22 deletions(-) create mode 100644 crates/resolved/tests/gpu/noise.rs diff --git a/crates/resolved-cli/src/main.rs b/crates/resolved-cli/src/main.rs index 54f6cc8..4d6b819 100644 --- a/crates/resolved-cli/src/main.rs +++ b/crates/resolved-cli/src/main.rs @@ -59,14 +59,21 @@ struct RenderArgs { lighting: LightingTier, /// Camera fidelity tier: C0, or C1 (adds pixels that collect light over their area, /// for graded edges and continuous thin parts, lens distortion, lens blur, and the - /// Sun, Moon and Earth where light time and aberration make them appear, and at - /// lighting L1 stars). + /// Sun, Moon and Earth where light time and aberration make them appear, detector + /// noise, and at lighting L1 stars). #[arg(long, default_value = "C0")] camera: CameraTier, /// From camera C1: render each pixel as N x N finer ones (1 to 4; default 2), which /// the lens blur then acts on before they are averaged into the pixel. #[arg(long, value_name = "N")] supersampling: Option, + /// From camera C1: the seed of the detector's noise. + #[arg(long, default_value_t = 0)] + seed: u64, + /// From camera C1: the frame number the noise is drawn for; the same seed and frame + /// give the same image. + #[arg(long, default_value_t = 0)] + frame: u64, /// Device: auto, gpu, cpu, or part of an adapter's name. #[arg(long, default_value = "auto")] adapter: String, @@ -221,6 +228,7 @@ fn render(args: &RenderArgs) -> Result<(), String> { dn_in_place: None, supersampling: args.supersampling, stars, + noise_seed: args.seed, }; let start = Instant::now(); let mut renderer = Renderer::new(camera, options).map_err(|e| e.to_string())?; @@ -236,6 +244,7 @@ fn render(args: &RenderArgs) -> Result<(), String> { let mut frame = None; for _ in 0..timings.unwrap_or(1).max(1) { let start = Instant::now(); + renderer.set_frame_number(args.frame); let rendered = renderer .submit(&scene, &observer) .and_then(|p| p.wait()) diff --git a/crates/resolved-cli/src/scene_file.rs b/crates/resolved-cli/src/scene_file.rs index f5ae6e6..d21d164 100644 --- a/crates/resolved-cli/src/scene_file.rs +++ b/crates/resolved-cli/src/scene_file.rs @@ -68,6 +68,12 @@ pub struct CameraFile { pub bias: Option, pub full_well: Option, pub adc_bits: Option, + /// From camera C1: read noise in electrons RMS. + pub read_noise: Option, + /// From camera C1: dark current in electrons per second. + pub dark_current: Option, + /// From camera C1: whether the electrons collected carry shot noise. + pub shot_noise: Option, /// The Gaussian blur's sigma in pixels, from camera C1: one for every channel, one /// per channel, or 0 for none. pub psf_sigma: Option, @@ -163,6 +169,11 @@ impl SceneFile { if let Some(bits) = c.adc_bits { camera.detector.adc_bits = bits; } + set(&mut camera.detector.read_noise, c.read_noise); + set(&mut camera.detector.dark_current, c.dark_current); + if let Some(shot) = c.shot_noise { + camera.detector.shot_noise = shot; + } if let Some([u, v]) = c.principal_point { camera.principal_point = Some(Vector2::new(u, v)); } @@ -441,13 +452,19 @@ mod tests { #[test] fn camera_overrides_apply() { let file: SceneFile = toml::from_str(&format!( - "{MINIMAL}\n[camera]\nwidth = 640\nexposure = 2e-3\nadc_bits = 10\n" + "{MINIMAL}\n[camera]\nwidth = 640\nexposure = 2e-3\nadc_bits = 10\nread_noise = 2.5\n\ + dark_current = 0\nshot_noise = false\n" )) .unwrap(); let camera = file.camera(); assert_eq!((camera.width, camera.height), (640, 1080)); assert_eq!(camera.exposure, 2e-3); assert_eq!(camera.detector.adc_bits, 10); + let d = camera.detector; + assert_eq!( + (d.read_noise, d.dark_current, d.shot_noise), + (2.5, 0.0, false) + ); for (line, psf) in [ ( diff --git a/crates/resolved/src/camera.rs b/crates/resolved/src/camera.rs index 1888190..8c7256b 100644 --- a/crates/resolved/src/camera.rs +++ b/crates/resolved/src/camera.rs @@ -21,6 +21,10 @@ pub struct Channel { } /// How the detector turns electrons into numbers. +/// +/// The noise fields apply from camera tier C1 up. A detector with no read noise, no +/// dark current and no shot noise ([`Detector::noiseless`]) gives noise-free images at +/// C1, for ground truth and labels. #[derive(Clone, Copy, Debug, PartialEq)] pub struct Detector { /// Electrons per DN. @@ -31,6 +35,13 @@ pub struct Detector { pub full_well: f64, /// Bits per pixel out of the ADC; output saturates at `2^adc_bits - 1`. pub adc_bits: u32, + /// Electrons RMS the readout adds to every pixel, every frame. + pub read_noise: f64, + /// Electrons per second each pixel makes on its own, at the detector's operating + /// temperature. They add to the light's, with their own shot noise. + pub dark_current: f64, + /// Whether the electrons collected vary as a Poisson count about their mean. + pub shot_noise: bool, } impl Detector { @@ -49,6 +60,34 @@ impl Detector { .round() .clamp(0.0, self.dn_max() as f64) as u16 } + + /// The same detector with no read noise, dark current or shot noise. + pub fn noiseless(self) -> Self { + Self { + read_noise: 0.0, + dark_current: 0.0, + shot_noise: false, + ..self + } + } + + /// Whether any noise or dark current is set. + pub fn is_noisy(&self) -> bool { + self.shot_noise || self.read_noise > 0.0 || self.dark_current > 0.0 + } + + /// The variance, in electrons squared, of one pixel's reading for a mean of + /// `electrons` from light over an exposure of `exposure` seconds, below the full + /// well: the shot noise of the light and of the dark current, the read noise, and + /// the ADC's rounding to whole DN (gain^2 / 12). + pub fn variance(&self, electrons: f64, exposure: f64) -> f64 { + let shot = if self.shot_noise { + electrons + self.dark_current * exposure + } else { + 0.0 + }; + shot + self.read_noise * self.read_noise + self.gain * self.gain / 12.0 + } } /// How the optics spread light from one point over the focal plane: the point spread @@ -223,6 +262,8 @@ impl CameraConfig { /// - Three Gaussian channels named `r`, `g`, `b`, centered at 620, 540 and 460 nm, /// 80 nm wide at half maximum, with 60 % peak QE. /// - 20,000 electron full well, 12-bit ADC, 5 electrons per DN, 64 DN bias. + /// - 5 electrons RMS read noise, 20 electrons per second dark current, and shot + /// noise: a CMOS sensor at room temperature. /// - 0.1 ms exposure: a white Lambertian surface facing the Sun fills about half the /// well (48-57 % across the channels), so Earth and the target are both unsaturated. /// - No lens distortion. @@ -255,6 +296,9 @@ impl CameraConfig { bias: 64.0, full_well: 20_000.0, adc_bits: 12, + read_noise: 5.0, + dark_current: 20.0, + shot_noise: true, }, near_clip: 0.01, mount: Pose::new(Vector3::zeros(), UnitQuaternion::identity()), @@ -306,6 +350,10 @@ impl CameraConfig { if !(1..=16).contains(&d.adc_bits) { return fail("ADC bits must be between 1 and 16"); } + let nonnegative = |x: f64| x.is_finite() && x >= 0.0; + if !nonnegative(d.read_noise) || !nonnegative(d.dark_current) { + return fail("read noise and dark current must be nonnegative"); + } Ok(()) } @@ -509,6 +557,21 @@ mod tests { assert_eq!(tight.dn(1e6), 4095); } + #[test] + fn detector_variance_adds_its_sources() { + // 1000 electrons of light and 20 e-/s of dark current over 0.5 s, read noise + // of 5 e-, gain 5 e-/DN: 1000 + 10 + 25 + 25/12. + let d = CameraConfig::example().detector; + assert_relative_eq!(d.variance(1000.0, 0.5), 1037.0 + 1.0 / 12.0, epsilon = 1e-9); + assert!(d.is_noisy()); + let quiet = d.noiseless(); + assert!(!quiet.is_noisy()); + assert_relative_eq!(quiet.variance(1000.0, 0.5), 25.0 / 12.0, epsilon = 1e-12); + let mut c = CameraConfig::example(); + c.detector.read_noise = -1.0; + assert!(c.validate().is_err()); + } + #[test] fn boresight_projects_to_principal_point() { let k = CameraConfig::example().intrinsics(); diff --git a/crates/resolved/src/fidelity.rs b/crates/resolved/src/fidelity.rs index 614d47c..6506490 100644 --- a/crates/resolved/src/fidelity.rs +++ b/crates/resolved/src/fidelity.rs @@ -50,7 +50,10 @@ pub enum CameraTier { /// ([`Scene::velocities`](crate::Scene::velocities)). Target meshes are covered at several /// points per rendered pixel, and rendered pixels on the edge of Earth, the Moon or /// the Sun cast a grid of rays, 4 x 4 and more for bodies only a few pixels across. - /// More of the camera joins this tier as it is modeled. + /// The detector adds its shot, dark and read noise + /// ([`Detector`](crate::camera::Detector)), seeded by + /// [`RendererOptions::noise_seed`](crate::RendererOptions::noise_seed) and the + /// frame number. More of the camera joins this tier as it is modeled. C1, } diff --git a/crates/resolved/src/gpu/sensor.rs b/crates/resolved/src/gpu/sensor.rs index 72ae88e..37abf35 100644 --- a/crates/resolved/src/gpu/sensor.rs +++ b/crates/resolved/src/gpu/sensor.rs @@ -1,8 +1,10 @@ //! The pass that turns channel radiance into the camera's DN. //! //! It is the last pass of a frame: a compute shader that reads the radiance image and -//! writes [`FrameTargets::dn`]. No noise yet; the result for each pixel is -//! [`Detector::dn`](crate::camera::Detector::dn) of [`Radiometry::electrons`]. +//! writes [`FrameTargets::dn`]. Without noise the result for each pixel is +//! [`Detector::dn`](crate::camera::Detector::dn) of [`Radiometry::electrons`]; with +//! it, the electrons vary as [`Detector::variance`](crate::camera::Detector::variance) +//! gives. use super::profile::PassTimestamps; use super::targets::FrameTargets; @@ -18,6 +20,8 @@ pub struct SensorUniforms { intrinsics: [f32; 4], detector: [f32; 4], dims: [u32; 4], + noise: [f32; 4], + key: [u32; 4], distortion: DistortionUniform, } @@ -43,9 +47,22 @@ impl SensorUniforms { camera.channels.len() as u32, d.dn_max(), ], + noise: [ + if d.shot_noise { 1.0 } else { 0.0 }, + d.read_noise, + d.dark_current * camera.exposure, + 0.0, + ] + .map(|x| x as f32), + key: [0; 4], distortion: DistortionUniform::new(k.distortion.as_ref().filter(|_| distorted)), } } + + /// Keys the noise to frame `frame` of the sequence `seed` starts. + pub fn set_key(&mut self, seed: u64, frame: u64) { + self.key = [frame as u32, seed as u32, (seed >> 32) as u32, 0]; + } } const WORKGROUP: u32 = 256; @@ -62,16 +79,16 @@ pub struct SensorPass { } impl SensorPass { - /// Compiles the shader, adding the stars' electrons if `stars`, and allocates the - /// uniform buffer. - pub fn new(ctx: &GpuContext, stars: bool) -> Self { + /// Compiles the shader, adding the stars' electrons if `stars` and the detector's + /// noise if `noise`, and allocates the uniform buffer. + pub fn new(ctx: &GpuContext, stars: bool, noise: bool) -> Self { let module = ctx .device .create_shader_module(wgpu::ShaderModuleDescriptor { label: Some("sensor"), source: wgpu::ShaderSource::Wgsl( [ - &format!("const STARS: bool = {stars};\n"), + &format!("const STARS: bool = {stars};\nconst NOISE: bool = {noise};\n"), include_str!("../shaders/distortion.wgsl"), include_str!("../shaders/sensor.wgsl"), ] diff --git a/crates/resolved/src/render.rs b/crates/resolved/src/render.rs index be6969e..bf8c4d7 100644 --- a/crates/resolved/src/render.rs +++ b/crates/resolved/src/render.rs @@ -78,6 +78,10 @@ pub struct RendererOptions { /// The stars, drawn at lighting tier L1 and up when given, by camera tier C1 and up /// (see [`StarCatalog::load`]); ignored below them. pub stars: Option>, + /// From camera C1, the seed of the detector's noise. A frame's noise depends only + /// on the seed and its number ([`Frame::number`]), so the same seed and number give + /// the same image, and different numbers independent noise. + pub noise_seed: u64, } /// One rendered image. @@ -96,6 +100,9 @@ pub struct Frame { /// GPU time per pass, if [`RendererOptions::profile`] was set and the adapter /// supports timestamps. pub timings: Option, + /// The frame's number: counted from 0 by each renderer, or set by + /// [`Renderer::set_frame_number`]. + pub number: u64, } impl Frame { @@ -135,6 +142,9 @@ pub struct Renderer { /// their pass. stars: Option<(Arc, StarPhotometry, StarPass)>, sensor_uniforms: SensorUniforms, + noise_seed: u64, + /// The next frame's number. + frame: u64, models: Vec, /// The buffer the DN image is read from: the sensor pass's own output, or a staging /// copy of it (see [`FrameTargets::dn_mappable`]). @@ -312,7 +322,7 @@ impl Renderer { blur, shadows: ShadowPass::new(&ctx, options.shadow_map_size.unwrap_or(4096)), earthshine: EarthshinePass::new(&ctx), - sensor: SensorPass::new(&ctx, stars.is_some()), + sensor: SensorPass::new(&ctx, stars.is_some(), area && camera.detector.is_noisy()), stars: stars.map(|catalog| { ( catalog, @@ -321,6 +331,8 @@ impl Renderer { ) }), sensor_uniforms: SensorUniforms::new(&camera, &radiometry, warp.is_some()), + noise_seed: options.noise_seed, + frame: 0, distorted: warp.is_some(), render_intrinsics, ctx, @@ -413,6 +425,12 @@ impl Renderer { Ok(()) } + /// Numbers the next frame submitted `number`, and those after it on from there: + /// to render a frame again with the same noise. + pub fn set_frame_number(&mut self, number: u64) { + self.frame = number; + } + /// The tiers this renderer draws at. pub fn fidelity(&self) -> Fidelity { self.fidelity @@ -593,6 +611,7 @@ impl Renderer { } None => empty(&mut encoder, "stars (none)", 5), } + self.sensor_uniforms.set_key(self.noise_seed, self.frame); self.sensor.encode( ctx, &mut encoder, @@ -641,10 +660,13 @@ impl Renderer { .slice(..) .map_async(wgpu::MapMode::Read, move |r| sender.send(r).unwrap_or(())); } + let number = self.frame; + self.frame += 1; Ok(Pending { renderer: Some(self), index: Some(index), mapped, + number, }) } @@ -672,6 +694,7 @@ pub struct Pending<'a> { renderer: Option<&'a mut Renderer>, index: Option, mapped: mpsc::Receiver>, + number: u64, } impl<'a> Pending<'a> { @@ -755,6 +778,7 @@ impl<'a> Pending<'a> { channels: r.camera.channels.len() as u32, radiance, timings, + number: self.number, }) } @@ -807,6 +831,8 @@ pub struct FrameView<'a> { pub radiance: Option>, /// GPU time per pass, as in [`Frame::timings`]. pub timings: Option, + /// The frame's number, as in [`Frame::number`]. + pub number: u64, } impl FrameView<'_> { @@ -831,6 +857,7 @@ impl FrameView<'_> { dn: self.dn().to_vec(), radiance: self.radiance.take(), timings: self.timings, + number: self.number, } } } diff --git a/crates/resolved/src/shaders/sensor.wgsl b/crates/resolved/src/shaders/sensor.wgsl index 3c59274..dc8e78e 100644 --- a/crates/resolved/src/shaders/sensor.wgsl +++ b/crates/resolved/src/shaders/sensor.wgsl @@ -1,5 +1,10 @@ // Channel radiance to DN: the camera equation, then the detector. // +// With NOISE (declared by the code that includes this file, from camera C1 with a +// noisy detector): shot noise of the light and dark current, then read noise, drawn +// from a counter-based generator keyed by the seed, the frame number and the value's +// index, so a frame's noise is the same every time it is rendered. +// // Output is u16 per channel, channels interleaved per pixel, rows from the top. WGSL // has no 16-bit integers, so each invocation computes two consecutive values and packs // them into one u32 word, low half first. Two invocations never touch the same word. @@ -11,6 +16,11 @@ struct Sensor { detector: vec4, // Width, height, channel count, largest DN. dims: vec4, + // Shot noise (1 or 0), read noise in electrons RMS, dark electrons over the + // exposure. + noise: vec4, + // The frame number, then the seed's low and high words. + key: vec4, // The lens distortion (distortion.wgsl, included first). distortion: Distortion, } @@ -22,6 +32,45 @@ struct Sensor { // value, 1/64 of an electron per unit (stars.wgsl). @group(0) @binding(3) var star_electrons: array; +// A hash of four words into four (Jarzynski and Olano 2020, pcg4d). +fn pcg4d(input: vec4) -> vec4 { + var v = input * 1664525u + 1013904223u; + v.x += v.y * v.w; + v.y += v.z * v.x; + v.z += v.x * v.y; + v.w += v.y * v.z; + v ^= v >> vec4(16u); + v.x += v.y * v.w; + v.y += v.z * v.x; + v.z += v.x * v.y; + v.w += v.y * v.z; + return v; +} + +// Uniform in (0, 1), from the top 24 bits. +fn unit(x: u32) -> f32 { + return (f32(x >> 8u) + 0.5) / 16777216.0; +} + +// A Poisson count of mean `mean` (at most SMALL) by inverting its distribution at `u`. +fn poisson(mean: f32, u: f32) -> f32 { + var p = exp(-mean); + var cdf = p; + var k = 0.0; + // Past the mean, the search ends once the terms no longer change the sum: for u + // within f32's rounding of 1, which the sum may never reach (one value in 196,608 + // at a mean of 3), it returns a count in the far tail. + while u > cdf && (k < mean || p > cdf * 6e-8) { + k += 1.0; + p *= mean / k; + cdf += p; + } + return k; +} + +// Above this mean a Poisson count is drawn as a Gaussian of the same variance. +const SMALL: f32 = 50.0; + fn value(i: u32) -> u32 { let n = s.dims.z; if i >= s.dims.x * s.dims.y * n { @@ -49,7 +98,23 @@ fn value(i: u32) -> u32 { // Light from the stars adds before the pixel fills. collected += f32(star_electrons[i]) / 64.0; } - let electrons = clamp(collected, 0.0, s.detector.y); + var electrons = clamp(collected, 0.0, s.detector.y); + if NOISE { + let r = pcg4d(vec4(i, s.key.x, s.key.y, s.key.z)); + // Two independent standard normals (Box-Muller). + let radius = sqrt(-2.0 * log(unit(r.x))); + let angle = 6.2831853 * unit(r.y); + let mean = collected + s.noise.z; + var count = mean; + if s.noise.x > 0.5 { + if mean > SMALL { + count = mean + sqrt(mean) * radius * cos(angle); + } else { + count = poisson(mean, unit(r.z)); + } + } + electrons = clamp(count, 0.0, s.detector.y) + s.noise.y * radius * sin(angle); + } let level = floor(electrons / s.detector.z + s.detector.w + 0.5); return u32(clamp(level, 0.0, f32(s.dims.w))); } diff --git a/crates/resolved/tests/gpu/distortion.rs b/crates/resolved/tests/gpu/distortion.rs index f943fb1..3119ade 100644 --- a/crates/resolved/tests/gpu/distortion.rs +++ b/crates/resolved/tests/gpu/distortion.rs @@ -302,7 +302,7 @@ fn the_sensor_uses_the_undistorted_line_of_sight() { targets.radiance.size(), ); let mut encoder = ctx.device.create_command_encoder(&Default::default()); - SensorPass::new(ctx, false).encode( + SensorPass::new(ctx, false, false).encode( ctx, &mut encoder, &targets, diff --git a/crates/resolved/tests/gpu/main.rs b/crates/resolved/tests/gpu/main.rs index e241924..dccb818 100644 --- a/crates/resolved/tests/gpu/main.rs +++ b/crates/resolved/tests/gpu/main.rs @@ -14,6 +14,7 @@ mod eclipses; mod gpu_context; mod mesh; mod moon_texture; +mod noise; mod pixel_area; mod psf; mod renderer; diff --git a/crates/resolved/tests/gpu/noise.rs b/crates/resolved/tests/gpu/noise.rs new file mode 100644 index 0000000..6d88886 --- /dev/null +++ b/crates/resolved/tests/gpu/noise.rs @@ -0,0 +1,364 @@ +//! Camera tier C1's detector noise against its statistics: the photon transfer curve, +//! dark frames, exact Poisson counts, reproducibility, and a faint star's SNR, on the +//! GPU and on the CPU (lavapipe). + +use crate::common; + +use std::sync::Arc; + +use resolved::camera::Detector; +use resolved::gpu::sensor::{SensorPass, SensorUniforms}; +use resolved::gpu::targets::{FrameTargets, TargetLayout}; +use resolved::gpu::{AdapterChoice, GpuContext}; +use resolved::nalgebra::{UnitQuaternion, Vector2, Vector3}; +use resolved::radiometry::{Radiometry, AU}; +use resolved::stars::CatalogStar; +use resolved::{ + CameraConfig, CameraTier, EarthTexture, Fidelity, LightingTier, Pose, Renderer, + RendererOptions, Scene, StarCatalog, StarPhotometry, +}; + +const SIZE: u32 = 256; + +/// 256 x 256 pixels of a 1 m focal length, so the falloff off axis is under 1e-5, with +/// `detector`. +fn camera(detector: Detector) -> CameraConfig { + CameraConfig { + width: SIZE, + height: SIZE, + focal_length: 1.0, + detector, + ..CameraConfig::example() + } +} + +/// The example detector with its noise replaced. +fn detector(read_noise: f64, dark_current: f64, shot_noise: bool) -> Detector { + Detector { + read_noise, + dark_current, + shot_noise, + ..CameraConfig::example().detector + } +} + +/// A detector that reads whole electrons: 1 e-/DN, 16 bits, a bias of 100 DN. +fn counting(read_noise: f64, dark_current: f64) -> Detector { + Detector { + gain: 1.0, + bias: 100.0, + full_well: 60_000.0, + adc_bits: 16, + read_noise, + dark_current, + shot_noise: true, + } +} + +/// Runs the sensor pass with noise over a flat field that gives `electrons` in every +/// pixel and channel, for frame `frame` of seed `seed`, and returns the DN image. +fn flat( + ctx: &GpuContext, + config: &CameraConfig, + electrons: f64, + seed: u64, + frame: u64, +) -> Vec { + let targets = FrameTargets::new(ctx, &TargetLayout::new(SIZE, SIZE, 3)); + let radiometry = Radiometry::new(config); + let level = (electrons / radiometry.electrons_per_radiance) as f32; + let radiance = vec![[level, level, level, 0.0]; (SIZE * SIZE) as usize]; + ctx.queue.write_texture( + wgpu::TexelCopyTextureInfo { + texture: &targets.radiance, + mip_level: 0, + origin: wgpu::Origin3d::ZERO, + aspect: wgpu::TextureAspect::All, + }, + bytemuck::cast_slice(&radiance), + wgpu::TexelCopyBufferLayout { + offset: 0, + bytes_per_row: Some(SIZE * 16), + rows_per_image: Some(SIZE), + }, + targets.radiance.size(), + ); + let pass = SensorPass::new(ctx, false, true); + let mut uniforms = SensorUniforms::new(config, &radiometry, false); + uniforms.set_key(seed, frame); + let mut encoder = ctx.device.create_command_encoder(&Default::default()); + pass.encode(ctx, &mut encoder, &targets, &uniforms, None); + ctx.queue.submit([encoder.finish()]); + targets.read_dn(ctx) +} + +fn mean(x: &[f64]) -> f64 { + x.iter().sum::() / x.len() as f64 +} + +fn variance(x: &[f64]) -> f64 { + let m = mean(x); + x.iter().map(|v| (v - m) * (v - m)).sum::() / (x.len() - 1) as f64 +} + +fn correlation(a: &[f64], b: &[f64]) -> f64 { + let (ma, mb) = (mean(a), mean(b)); + let cov: f64 = a.iter().zip(b).map(|(x, y)| (x - ma) * (y - mb)).sum(); + cov / (variance(a) * variance(b)).sqrt() / (a.len() - 1) as f64 +} + +/// Channel `k`'s values of a DN image, as f64. +fn channel(dn: &[u16], k: usize) -> Vec { + dn.iter().skip(k).step_by(3).map(|&d| d as f64).collect() +} + +#[test] +fn flat_fields_follow_the_photon_transfer_curve() { + // Two frames of a flat field at each of six levels, as EMVA 1288 measures a camera: + // the mean of the two less the bias, and half the variance of their difference, + // which leaves out anything the same in both. The temporal variance in DN^2 is + // (electrons + read noise^2) / gain^2 + 1/12 for the rounding, and the GPU and + // lavapipe each match it to 3 % (65,536 pixels give a sampling error of 0.6 %), + // and a line through the levels has slope 1 / gain. + let config = camera(detector(5.0, 0.0, true)); + let d = config.detector; + for ctx in common::contexts() { + let name = &ctx.info.name; + let mut points = Vec::new(); + for (level, electrons) in [100.0, 500.0, 2000.0, 5000.0, 10_000.0, 16_000.0] + .into_iter() + .enumerate() + { + // Frames of their own at each level, so the levels sample independently. + let a = flat(ctx, &config, electrons, 7, 2 * level as u64); + let b = flat(ctx, &config, electrons, 7, 2 * level as u64 + 1); + for k in 0..3 { + let (a, b) = (channel(&a, k), channel(&b, k)); + let m = mean( + &a.iter() + .zip(&b) + .map(|(x, y)| (x + y) / 2.0) + .collect::>(), + ) - d.bias; + let v = variance(&a.iter().zip(&b).map(|(x, y)| x - y).collect::>()) / 2.0; + let want = d.variance(electrons, config.exposure) / (d.gain * d.gain); + assert!( + (m / (electrons / d.gain) - 1.0).abs() < 2e-3, + "{name}: {electrons} e-, channel {k}: mean {m} DN" + ); + assert!( + (v / want - 1.0).abs() < 0.03, + "{name}: {electrons} e-, channel {k}: variance {v} DN^2, expected {want}" + ); + points.push((m, v)); + } + } + let (mx, my) = ( + mean(&points.iter().map(|p| p.0).collect::>()), + mean(&points.iter().map(|p| p.1).collect::>()), + ); + let slope = points.iter().map(|p| (p.0 - mx) * (p.1 - my)).sum::() + / points.iter().map(|p| (p.0 - mx) * (p.0 - mx)).sum::(); + assert!( + (slope * d.gain - 1.0).abs() < 0.02, + "{name}: slope {slope} DN, expected {}", + 1.0 / d.gain + ); + } +} + +#[test] +fn dark_frames_carry_dark_current_and_read_noise() { + // No light, 1000 e-/s of dark current over 0.1 s: a mean of 100 electrons, 20 DN + // above the bias, and a variance of (100 + 5^2) / 25 + 1/12 DN^2. Without shot + // noise the dark electrons are only an offset, and the variance is the read + // noise's and the rounding's. + for ctx in common::contexts() { + let name = &ctx.info.name; + for shot in [true, false] { + let config = CameraConfig { + exposure: 0.1, + ..camera(detector(5.0, 1000.0, shot)) + }; + let d = config.detector; + let dn = channel(&flat(ctx, &config, 0.0, 3, 0), 1); + let m = mean(&dn) - d.bias; + let v = variance(&dn); + let want = d.variance(0.0, config.exposure) / (d.gain * d.gain); + assert!((m - 20.0).abs() < 0.05, "{name}: shot {shot}: mean {m} DN"); + assert!( + (v / want - 1.0).abs() < 0.03, + "{name}: shot {shot}: variance {v} DN^2, expected {want}" + ); + } + } +} + +#[test] +fn small_counts_are_poisson() { + // At 1 e-/DN with no read noise the DN less the bias is the count itself. Below a + // mean of 50 it is drawn exactly: every count's frequency over 196,608 values + // matches the Poisson probability to 5 standard errors. Above 50 the Gaussian of + // the same variance takes over: its mean and variance still match. + for ctx in common::contexts() { + let name = &ctx.info.name; + let config = camera(counting(0.0, 0.0)); + for mean_count in [0.5, 3.0, 20.0, 45.0, 60.0] { + let dn = flat(ctx, &config, mean_count, 11, 0); + let counts: Vec = dn.iter().map(|&d| d as f64 - 100.0).collect(); + let n = counts.len() as f64; + let (m, v) = (mean(&counts), variance(&counts)); + assert!( + (m / mean_count - 1.0).abs() < 0.02 && (v / mean_count - 1.0).abs() < 0.03, + "{name}: mean {mean_count}: measured mean {m}, variance {v}" + ); + if mean_count > 50.0 { + continue; + } + let mut p = (-mean_count).exp(); + for k in 0..(mean_count * 3.0 + 10.0) as usize { + if k > 0 { + p *= mean_count / k as f64; + } + let seen = counts.iter().filter(|&&c| c == k as f64).count() as f64 / n; + let error = (p * (1.0 - p) / n).sqrt(); + assert!( + (seen - p).abs() < 5.0 * error + 1e-5, + "{name}: mean {mean_count}: P({k}) {seen}, Poisson {p}" + ); + } + } + } +} + +#[test] +fn noise_repeats_by_seed_and_frame_and_is_independent() { + // The same seed and frame give the same image; another frame or another seed gives + // noise uncorrelated with it (|r| < 0.01 over 196,608 values, against a sampling + // error of 0.0023), and neighboring values within a frame are uncorrelated too. + let config = camera(detector(5.0, 0.0, true)); + for ctx in common::contexts() { + let name = &ctx.info.name; + let base = flat(ctx, &config, 1000.0, 1, 0); + assert_eq!(base, flat(ctx, &config, 1000.0, 1, 0), "{name}: repeat"); + let values: Vec = base.iter().map(|&d| d as f64).collect(); + for (what, other) in [ + ("next frame", flat(ctx, &config, 1000.0, 1, 1)), + ("next seed", flat(ctx, &config, 1000.0, 2, 0)), + ( + "seed's high word", + flat(ctx, &config, 1000.0, 1 + (1 << 32), 0), + ), + ] { + let other: Vec = other.iter().map(|&d| d as f64).collect(); + let r = correlation(&values, &other); + assert!(r.abs() < 0.01, "{name}: {what}: r = {r}"); + } + let r = correlation(&values[..values.len() - 1], &values[1..]); + assert!(r.abs() < 0.01, "{name}: neighboring values: r = {r}"); + let r = correlation(&values[..values.len() - 3], &values[3..]); + assert!(r.abs() < 0.01, "{name}: neighboring pixels: r = {r}"); + } +} + +#[test] +fn a_faint_star_has_the_signal_to_noise_the_detector_gives() { + // 400 identical stars of V 11.5 in a 20 x 20 grid, 12 pixels apart, through the + // example detector's noise, exposed to give each 300 electrons. Each star's electrons summed over the 7 x 7 + // pixels around it, less the bias, average S + 49 I_d t, S its electrons and I_d t + // each pixel's dark electrons, and vary from star to star with variance + // S + 49 (I_d t + read^2 + gain^2 / 12) electrons^2: an SNR of about 7. The sums' + // mean and spread match that to 2 % and 12 % (400 stars measure a spread to 3.5 %). The + // renderer numbers its frames, and setting the number back repeats a frame. + let _devices = common::device_lock(); + let photometry = StarPhotometry::new(&camera(CameraConfig::example().detector)); + let signal = 300.0; + let exposure = signal / (photometry.rate(1, 11.5, 0.6) * photometry.aperture); + let config = CameraConfig { + exposure, + ..camera(CameraConfig::example().detector) + }; + let k = config.intrinsics(); + let position = Vector3::new(4.2e7, 0.0, 0.0); + let look = Pose::looking_at(position, position + Vector3::y(), Vector3::z()).unwrap(); + let mut stars = Vec::new(); + let mut centers = Vec::new(); + for j in 0..20 { + for i in 0..20 { + let uv = Vector2::new(14.3 + 12.0 * i as f64, 14.7 + 12.0 * j as f64); + let d = look.attitude * k.ray(&uv).normalize(); + stars.push(CatalogStar { + ra: d.y.atan2(d.x), + dec: d.z.asin(), + pm_ra: 0.0, + pm_dec: 0.0, + v: 11.5, + bv: 0.6, + }); + centers.push(uv); + } + } + let d = config.detector; + let dark = d.dark_current * exposure; + let per_pixel = d.variance(0.0, exposure); + let spread = (signal + 49.0 * per_pixel).sqrt(); + let scene = Scene { + earth_attitude: UnitQuaternion::identity(), + sun: position - Vector3::y() * AU, + moon: Vector3::new(0.0, 0.0, -3.844e8), + moon_attitude: UnitQuaternion::identity(), + targets: Vec::new(), + velocities: Default::default(), + }; + let catalog = Arc::new(StarCatalog::new(stars, 0.0)); + for adapter in [AdapterChoice::Gpu, AdapterChoice::Cpu] { + let black = EarthTexture::from_rgba8(1, 1, vec![0, 0, 0, 255]).unwrap(); + let options = RendererOptions { + adapter, + fidelity: Fidelity { + lighting: LightingTier::L1, + camera: CameraTier::C1, + }, + earth_texture: Some(Arc::new(black)), + moon_texture: Some(common::uniform_moon()), + stars: Some(catalog.clone()), + noise_seed: 5, + ..Default::default() + }; + let mut r = Renderer::new(config.clone(), options).unwrap(); + let name = r.adapter_info().name.clone(); + let first = r.submit(&scene, &look).unwrap().wait().unwrap(); + let second = r.submit(&scene, &look).unwrap().wait().unwrap(); + assert_eq!((first.number, second.number), (0, 1), "{name}"); + assert_ne!(first.dn, second.dn, "{name}"); + r.set_frame_number(0); + let again = r.submit(&scene, &look).unwrap().wait().unwrap(); + assert_eq!(again.number, 0, "{name}"); + assert_eq!(first.dn, again.dn, "{name}: frame 0 again"); + + let sums: Vec = centers + .iter() + .map(|c| { + let (u0, v0) = (c.x as u32, c.y as u32); + let mut sum = 0.0; + for v in v0 - 3..=v0 + 3 { + for u in u0 - 3..=u0 + 3 { + sum += (first.dn_at(u, v, 1) as f64 - d.bias) * d.gain; + } + } + sum + }) + .collect(); + let (m, s) = (mean(&sums), variance(&sums).sqrt()); + let want = signal + 49.0 * dark; + assert!( + (m / want - 1.0).abs() < 0.02, + "{name}: mean {m} e-, expected {want}" + ); + assert!( + (s / spread - 1.0).abs() < 0.12, + "{name}: spread {s} e-, expected {spread} (SNR {})", + signal / spread + ); + } +} diff --git a/crates/resolved/tests/gpu/sensor.rs b/crates/resolved/tests/gpu/sensor.rs index 5f0a5d9..eec3451 100644 --- a/crates/resolved/tests/gpu/sensor.rs +++ b/crates/resolved/tests/gpu/sensor.rs @@ -32,7 +32,7 @@ fn sense(ctx: &GpuContext, name: &str, radiance: &[[f32; 4]]) -> Vec { }, targets.radiance.size(), ); - let pass = SensorPass::new(ctx, false); + let pass = SensorPass::new(ctx, false, false); let uniforms = SensorUniforms::new(&config, &Radiometry::new(&config), false); let mut encoder = ctx.device.create_command_encoder(&Default::default()); pass.encode(ctx, &mut encoder, &targets, &uniforms, None); @@ -164,7 +164,7 @@ fn sunlit_earth_through_the_whole_chain() { let targets = FrameTargets::new(ctx, &TargetLayout::new(config.width, config.height, 3)); let mut encoder = ctx.device.create_command_encoder(&Default::default()); BodiesPass::new(ctx, 1, false).encode(ctx, &mut encoder, &targets, &bodies, None); - SensorPass::new(ctx, false).encode(ctx, &mut encoder, &targets, &sensor, None); + SensorPass::new(ctx, false, false).encode(ctx, &mut encoder, &targets, &sensor, None); ctx.queue.submit([encoder.finish()]); let dn = targets.read_dn(ctx); common::save_dn("sensor-earth", &ctx.info, &dn); diff --git a/crates/resolved/tests/gpu/stars.rs b/crates/resolved/tests/gpu/stars.rs index ab4e1ee..84dd848 100644 --- a/crates/resolved/tests/gpu/stars.rs +++ b/crates/resolved/tests/gpu/stars.rs @@ -36,6 +36,9 @@ fn camera(exposure: f64) -> CameraConfig { bias: 100.0, full_well: 60_000.0, adc_bits: 16, + read_noise: 0.0, + dark_current: 0.0, + shot_noise: false, }, ..CameraConfig::example() } diff --git a/docs/optics.md b/docs/optics.md index d195a4e..d85ee0f 100644 --- a/docs/optics.md +++ b/docs/optics.md @@ -356,7 +356,8 @@ For the example camera, pi * A * t / (4 N^2) is 1.48e-16 m^2 sr s. ## From electrons to DN: the detector -Code: `camera.rs` (`Detector::dn`); `gpu/sensor.rs`, `shaders/sensor.wgsl`. +Code: `camera.rs` (`Detector::dn`, `Detector::variance`); `gpu/sensor.rs`, +`shaders/sensor.wgsl`. - **Full well.** The maximum number of electrons a pixel holds. Beyond it the pixel saturates and reads its maximum value whatever the light. @@ -371,7 +372,8 @@ Code: `camera.rs` (`Detector::dn`); `gpu/sensor.rs`, `shaders/sensor.wgsl`. So DN = bias + min(electrons, full well) / gain, rounded and clipped to the ADC range. The renderer does this on the GPU for every pixel and channel, as the last step of a frame, and reads back only the resulting integers. -Noise, not modeled yet, enters between the electrons and the gain. +From camera tier C1, noise enters between the electrons and the gain (see +Noise below). Two consequences of these numbers show up in every image: @@ -411,7 +413,98 @@ exactly, on the GPU and on lavapipe. Picking the exposure time is the same trade a photographer makes. A longer exposure saturates the sunlit target. A shorter one leaves the side of the -target lit only by Earthshine lost in noise once noise is modeled. +target lit only by Earthshine lost in noise. + +### Noise + +Code: `Detector::variance`; `sensor.wgsl` (`NOISE`). + +A real image is noisy, and the noise decides what can be detected at all. +From camera tier C1 the renderer adds three kinds, following the EMVA 1288 +standard's model of a camera: + +- **Shot noise.** Light arrives as photons at random times, so a pixel + expecting N electrons collects a Poisson-distributed count with mean N and + standard deviation sqrt(N). It is a property of light, not of the camera: + a perfect detector still has it. A pixel expecting 10,000 electrons varies + by 100 (1 %); one expecting 100 varies by 10 (10 %). +- **Dark current.** The detector makes electrons of its own from heat, + `Detector::dark_current` per second, which add to the light's and carry + shot noise of their own. It grows with temperature, roughly doubling every + 6 to 8 degrees C; the value given is for the detector's operating + temperature. +- **Read noise.** Reading a pixel out adds a few electrons of random error, + `Detector::read_noise` RMS, to every pixel in every frame, light or not. + +In order: the light's electrons and the dark electrons together are drawn +as one Poisson count (a sum of Poisson counts is itself one); the count is +clipped to the full well; read noise is added; then the gain, bias and +rounding to DN as before. Below a mean of 50 electrons the count is drawn +exactly. Above it a Gaussian of the same mean and variance takes its place, +which a Poisson count's shape is already close to and which costs the same +at any level. With read noise the DN can fall below the bias, which is what +the bias is for. + +Each pixel's variance in electrons^2 is then + + variance = (electrons + dark current x exposure) + read noise^2 + gain^2 / 12 + +where the last term is the rounding to whole DN. Divided by gain^2 it is in +DN^2. Plotted against the mean signal in DN, the variance is a straight line +with slope 1 / gain and intercept (read noise^2 + dark shot noise) / gain^2: +the photon transfer curve, the standard way to measure a camera's gain and +read noise from flat fields alone. + +For the example camera (5 electrons read noise, 20 electrons per second dark +current, gain 5) at 0.1 ms, the signal-to-noise ratio, signal over the +square root of the variance, is: + +| What | Electrons | SNR | +|---|---|---| +| A white surface facing the Sun, green channel | 11,139 | 105 | +| Earth, albedo 0.3 | 3,342 | 58 | +| A white surface lit by 1 % of sunlight, as Earthshine can | 111 | 9.5 | +| A V = 0 star, over the 4 pixels it covers | 56 | 4.4 | +| A V = 4 star, over 4 pixels | 1.4 | 0.13 | + +So at the exposures a sunlit target needs, the brightest stars are just +detectable and the rest are lost in read noise. At 50 ms, a star tracker's +exposure, a V = 6 star reaches an SNR of 7.5. + +Dark current hardly matters at these exposures: 20 electrons per second is +0.002 electrons in 0.1 ms. It matters for long ones: 1 electron in 50 ms, +100 in 5 s, with 10 electrons of shot noise. + +**Reproducible randomness.** The noise comes from a counter-based generator: +a hash (Jarzynski and Olano's pcg4d) of the pixel and channel, the frame +number and a 64-bit seed (`RendererOptions::noise_seed`), turned into +uniform numbers, then into Gaussians (Box-Muller) or a Poisson count (by +inverting its distribution). Nothing is carried from one pixel or frame to +the next, so every value is computed independently on the GPU. The renderer +numbers its frames from 0 (`Frame::number`), and +`Renderer::set_frame_number` renders a frame again with exactly the same +noise. Different frame numbers or seeds give independent noise. The hash is +identical on every device, but the logarithms and cosines that shape it may +differ in the last bit, so the GPU and lavapipe agree statistically rather +than pixel for pixel. + +**Noise-free images.** A detector with no read noise, no dark current and +shot noise off (`Detector::noiseless`) renders C1 without noise, exactly as +before noise existed, for ground truth and ML labels. The renderer then +compiles the noise out entirely. At C0 there is never noise. + +Measured, identically on the GPU and lavapipe: flat fields at six levels +from 100 to 16,000 electrons, two frames each, give half the variance of +their difference within 2 % of the formula (0.982 to 1.008 of it, against a +sampling error of 0.6 %), and the photon transfer curve's slope is 1 / gain +to 0.1 %. A dark frame of 100 electrons reads 20.005 DN above the bias, with +its variance within 0.6 %. Counts below 50 match the Poisson probabilities +count by count. Consecutive frames, different seeds and neighboring pixels +are uncorrelated (|r| at most 0.005 over 196,608 values). 400 faint stars of +300 electrons each vary as predicted to 0.4 %, an SNR of 7.35. Noise costs nothing measurable on an +RTX 3060 (0.476 ms for the sensor pass at 1920 x 1080 with or without it, +1.89 ms at 3840 x 2160); on lavapipe it adds 1.5 ms at 1920 x 1080 and +9.6 ms at 3840 x 2160. ## How shapes become pixels @@ -1461,7 +1554,7 @@ is a run of the script with `--gaia 16`, at twice the size. - Earthshine on the Moon's night side - The diffraction pattern's rings, and a lens's measured point spread function -- Noise: photon shot noise, read noise, dark current, pixel-to-pixel variation +- Pixel-to-pixel variation in response and dark current (PRNU, DSNU) - Saturation spilling into neighbors (blooming) - Non-Lambertian materials (BRDFs), glints, and the full Moon's edge as bright as its center diff --git a/readme.md b/readme.md index da6524c..8643301 100644 --- a/readme.md +++ b/readme.md @@ -173,6 +173,8 @@ cargo run --release -p resolved-showcase -- --serve 0.0.0.0:8010 L0 and C0). `--supersampling N` renders each pixel as N x N finer ones from C1 (default 2). A tier the scene or the data cannot support stops with an error naming what is missing. +- `--seed N` and `--frame N` pick the detector noise at camera C1 (default 0 + and 0): the same seed and frame give the same image. - `--data DIR` says where the tiers' data is (`earth/bmng-2004-MM.jpg`, `moon/wac-hapke-643.png` and `stars/stars.bin`). Without it the CLI tries `RESOLVED_DATA`, then `./data`, then the @@ -202,10 +204,11 @@ Earth and the Moon, the Sun's limb darkening, and stars from Hipparcos, Tycho-2 and Gaia, which camera C1 shows. Camera C1 is under way: so far, pixels that collect light over their area, lens distortion (Brown-Conrady), lens blur (a Gaussian PSF) applied to the scene rendered 2 x 2 finer than the -pixels, and the Sun, Moon and Earth at their apparent positions (light time -and aberration). Planned, in order: +pixels, the Sun, Moon and Earth at their apparent positions (light time +and aberration), and detector noise (shot, dark current and read noise, +seeded and reproducible). Planned, in order: -- Camera: detector noise, pixel-to-pixel variation, saturation and blooming. +- Camera: pixel-to-pixel variation, saturation and blooming. - Lighting: BRDF and spectral materials, atmospheric scattering, clouds, lunar terrain. - Portable ray tracing, path tracing, and label outputs for ML.