Sensor/detector DNs

crates/resolved/src/camera.rs

@@ -33,6 +33,24 @@ pub struct Detector {
     pub adc_bits: u32,
 }
 
+impl Detector {
+    /// The largest value the ADC can output.
+    pub fn dn_max(&self) -> u32 {
+        (1 << self.adc_bits) - 1
+    }
+
+    /// The DN a pixel reads after collecting `electrons`, with no noise: the charge
+    /// saturates at the full well, is divided by the gain, offset by the bias, rounded
+    /// to the nearest integer and clipped to the ADC's range. The sensor pass on the GPU
+    /// does the same in f32.
+    pub fn dn(&self, electrons: f64) -> u16 {
+        let charge = electrons.clamp(0.0, self.full_well);
+        (charge / self.gain + self.bias)
+            .round()
+            .clamp(0.0, self.dn_max() as f64) as u16
+    }
+}
+
 /// A pinhole camera with radiometric optics and a detector.
 #[derive(Clone, Debug, PartialEq)]
 pub struct CameraConfig {
@@ -232,6 +250,19 @@ mod tests {
         assert!(camera.validate().is_err());
     }
 
+    #[test]
+    fn detector_saturates_at_full_well_then_adc() {
+        let d = CameraConfig::example().detector;
+        assert_eq!(d.dn(0.0), 64);
+        assert_eq!(d.dn(-5.0), 64);
+        assert_eq!(d.dn(11_138.8), 2292);
+        // The full well fills before the ADC tops out: 20,000 / 5 + 64 = 4064 < 4095.
+        assert_eq!(d.dn(20_000.0), 4064);
+        assert_eq!(d.dn(5e8), 4064);
+        let tight = Detector { gain: 1.0, ..d };
+        assert_eq!(tight.dn(1e6), 4095);
+    }
+
     #[test]
     fn boresight_projects_to_principal_point() {
         let k = CameraConfig::example().intrinsics();

crates/resolved/src/gpu/mod.rs

@@ -3,6 +3,7 @@
 pub mod bodies;
 pub mod context;
 pub mod mesh;
+pub mod sensor;
 pub mod targets;
 
 pub use context::{adapters, AdapterChoice, GpuContext};

crates/resolved/src/gpu/sensor.rs

@@ -0,0 +1,124 @@
+//! 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 target and
+//! writes [`FrameTargets::dn`]. No noise yet (M3); the result for each pixel is
+//! [`Detector::dn`](crate::camera::Detector::dn) of [`Radiometry::electrons`].
+
+use super::targets::FrameTargets;
+use super::GpuContext;
+use crate::camera::CameraConfig;
+use crate::radiometry::Radiometry;
+
+/// The shader's uniform block; the layout matches `Sensor` in `sensor.wgsl`.
+#[repr(C)]
+#[derive(Clone, Copy, Debug, PartialEq, bytemuck::Pod, bytemuck::Zeroable)]
+pub struct SensorUniforms {
+    intrinsics: [f32; 4],
+    detector: [f32; 4],
+    dims: [u32; 4],
+}
+
+impl SensorUniforms {
+    /// The constants for `camera`, whose channels `radiometry` was built from. They do
+    /// not change from frame to frame.
+    pub fn new(camera: &CameraConfig, radiometry: &Radiometry) -> Self {
+        let k = camera.intrinsics();
+        let d = &camera.detector;
+        Self {
+            intrinsics: [k.fx, k.fy, k.cx, k.cy].map(|x| x as f32),
+            detector: [
+                radiometry.electrons_per_radiance,
+                d.full_well,
+                d.gain,
+                d.bias,
+            ]
+            .map(|x| x as f32),
+            dims: [
+                camera.width,
+                camera.height,
+                camera.channels.len() as u32,
+                d.dn_max(),
+            ],
+        }
+    }
+}
+
+const WORKGROUP: u32 = 256;
+
+/// The compiled sensor pipeline and its uniform buffer.
+#[derive(Debug)]
+pub struct SensorPass {
+    pipeline: wgpu::ComputePipeline,
+    uniforms: wgpu::Buffer,
+}
+
+impl SensorPass {
+    /// Compiles the shader and allocates the uniform buffer.
+    pub fn new(ctx: &GpuContext) -> Self {
+        let module = ctx
+            .device
+            .create_shader_module(wgpu::ShaderModuleDescriptor {
+                label: Some("sensor"),
+                source: wgpu::ShaderSource::Wgsl(include_str!("../shaders/sensor.wgsl").into()),
+            });
+        let pipeline = ctx
+            .device
+            .create_compute_pipeline(&wgpu::ComputePipelineDescriptor {
+                label: Some("sensor"),
+                layout: None,
+                module: &module,
+                entry_point: Some("main"),
+                compilation_options: Default::default(),
+                cache: None,
+            });
+        let uniforms = ctx.device.create_buffer(&wgpu::BufferDescriptor {
+            label: Some("sensor uniforms"),
+            size: std::mem::size_of::<SensorUniforms>() as u64,
+            usage: wgpu::BufferUsages::UNIFORM | wgpu::BufferUsages::COPY_DST,
+            mapped_at_creation: false,
+        });
+        Self { pipeline, uniforms }
+    }
+
+    /// Uploads `uniforms` and records the pass into `encoder`, reading `targets`'
+    /// radiance and writing its DN buffer.
+    pub fn encode(
+        &self,
+        ctx: &GpuContext,
+        encoder: &mut wgpu::CommandEncoder,
+        targets: &FrameTargets,
+        uniforms: &SensorUniforms,
+    ) {
+        ctx.queue
+            .write_buffer(&self.uniforms, 0, bytemuck::bytes_of(uniforms));
+        let radiance = targets.radiance.create_view(&Default::default());
+        let bind_group = ctx.device.create_bind_group(&wgpu::BindGroupDescriptor {
+            label: Some("sensor"),
+            layout: &self.pipeline.get_bind_group_layout(0),
+            entries: &[
+                wgpu::BindGroupEntry {
+                    binding: 0,
+                    resource: wgpu::BindingResource::TextureView(&radiance),
+                },
+                wgpu::BindGroupEntry {
+                    binding: 1,
+                    resource: targets.dn.as_entire_binding(),
+                },
+                wgpu::BindGroupEntry {
+                    binding: 2,
+                    resource: self.uniforms.as_entire_binding(),
+                },
+            ],
+        });
+        // One invocation per u32 word, in rows of at most 65,535 workgroups.
+        let groups = (targets.dn.size() / 4).div_ceil(WORKGROUP as u64) as u32;
+        let (x, y) = (groups.min(65_535), groups.div_ceil(65_535));
+        let mut pass = encoder.begin_compute_pass(&wgpu::ComputePassDescriptor {
+            label: Some("sensor"),
+            timestamp_writes: None,
+        });
+        pass.set_pipeline(&self.pipeline);
+        pass.set_bind_group(0, &bind_group, &[]);
+        pass.dispatch_workgroups(x, y, 1);
+    }
+}

crates/resolved/src/gpu/targets.rs

@@ -21,11 +21,16 @@ pub struct FrameTargets {
     pub radiance: wgpu::Texture,
     /// Depth, [`DEPTH_FORMAT`].
     pub depth: wgpu::Texture,
+    /// Number of camera channels.
+    pub channels: u32,
+    /// The sensor pass's output: u16 DN, channels interleaved per pixel, rows from the
+    /// top, packed two to a u32 word.
+    pub dn: wgpu::Buffer,
 }
 
 impl FrameTargets {
-    /// Allocates targets of `width` x `height` pixels.
+    /// Allocates targets of `width` x `height` pixels with `channels` channels.
-    pub fn new(ctx: &GpuContext, width: u32, height: u32) -> Self {
+    pub fn new(ctx: &GpuContext, width: u32, height: u32, channels: u32) -> Self {
         let texture = |label, format, usage| {
             ctx.device.create_texture(&wgpu::TextureDescriptor {
                 label: Some(label),
@@ -45,14 +50,57 @@ impl FrameTargets {
         let usage = wgpu::TextureUsages::RENDER_ATTACHMENT
             | wgpu::TextureUsages::TEXTURE_BINDING
             | wgpu::TextureUsages::COPY_SRC;
+        let values = width as u64 * height as u64 * channels as u64;
         Self {
             width,
             height,
-            radiance: texture("radiance", RADIANCE_FORMAT, usage),
+            // Writable from the CPU so tests can supply radiance directly.
+            radiance: texture(
+                "radiance",
+                RADIANCE_FORMAT,
+                usage | wgpu::TextureUsages::COPY_DST,
+            ),
             depth: texture("depth", DEPTH_FORMAT, usage),
+            channels,
+            dn: ctx.device.create_buffer(&wgpu::BufferDescriptor {
+                label: Some("dn"),
+                size: values.div_ceil(2) * 4,
+                usage: wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_SRC,
+                mapped_at_creation: false,
+            }),
         }
     }
 
+    /// Number of u16 values in the DN image.
+    pub fn dn_len(&self) -> usize {
+        self.width as usize * self.height as usize * self.channels as usize
+    }
+
+    /// Reads the DN image back. Slow, like [`FrameTargets::read_radiance`]; the per-frame
+    /// readback is the renderer's.
+    pub fn read_dn(&self, ctx: &GpuContext) -> Vec<u16> {
+        let staging = ctx.device.create_buffer(&wgpu::BufferDescriptor {
+            label: Some("dn readback"),
+            size: self.dn.size(),
+            usage: wgpu::BufferUsages::MAP_READ | wgpu::BufferUsages::COPY_DST,
+            mapped_at_creation: false,
+        });
+        let mut encoder = ctx.device.create_command_encoder(&Default::default());
+        encoder.copy_buffer_to_buffer(&self.dn, 0, &staging, 0, self.dn.size());
+        ctx.queue.submit([encoder.finish()]);
+        staging
+            .slice(..)
+            .map_async(wgpu::MapMode::Read, |r| r.expect("readback map"));
+        ctx.device
+            .poll(wgpu::PollType::wait_indefinitely())
+            .expect("readback poll");
+        let view = staging
+            .slice(..)
+            .get_mapped_range()
+            .expect("readback range");
+        bytemuck::cast_slice::<u8, u16>(&view)[..self.dn_len()].to_vec()
+    }
+
     /// Reads the radiance target back, row by row from the top. Slow; for tests and
     /// debugging, not the per-frame path.
     pub fn read_radiance(&self, ctx: &GpuContext) -> Vec<[f32; 4]> {

crates/resolved/src/shaders/sensor.wgsl

@@ -0,0 +1,47 @@
+// Channel radiance to DN: the camera equation, then the detector.
+//
+// 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.
+
+struct Sensor {
+    // fx, fy, cx, cy in pixels.
+    intrinsics: vec4<f32>,
+    // Electrons per unit channel radiance on axis, full well, gain (e-/DN), bias (DN).
+    detector: vec4<f32>,
+    // Width, height, channel count, largest DN.
+    dims: vec4<u32>,
+}
+
+@group(0) @binding(0) var radiance: texture_2d<f32>;
+@group(0) @binding(1) var<storage, read_write> dn: array<u32>;
+@group(0) @binding(2) var<uniform> s: Sensor;
+
+fn value(i: u32) -> u32 {
+    let n = s.dims.z;
+    if i >= s.dims.x * s.dims.y * n {
+        return 0u;
+    }
+    let pixel = i / n;
+    let x = pixel % s.dims.x;
+    let y = pixel / s.dims.x;
+    let l = textureLoad(radiance, vec2<u32>(x, y), 0)[i % n];
+
+    // Off-axis falloff: cos^4 of the angle between this pixel's line of sight and the
+    // boresight.
+    let k = s.intrinsics;
+    let ray = vec3<f32>((f32(x) + 0.5 - k.z) / k.x, (f32(y) + 0.5 - k.w) / k.y, 1.0);
+    let c2 = 1.0 / dot(ray, ray);
+    let electrons = clamp(l * s.detector.x * c2 * c2, 0.0, s.detector.y);
+    let level = floor(electrons / s.detector.z + s.detector.w + 0.5);
+    return u32(clamp(level, 0.0, f32(s.dims.w)));
+}
+
+@compute @workgroup_size(256)
+fn main(@builtin(global_invocation_id) id: vec3<u32>, @builtin(num_workgroups) groups: vec3<u32>) {
+    let word = id.x + id.y * groups.x * 256u;
+    if word >= arrayLength(&dn) {
+        return;
+    }
+    dn[word] = value(2u * word) | (value(2u * word + 1u) << 16u);
+}

crates/resolved/tests/bodies.rs

@@ -7,7 +7,7 @@ use approx::assert_relative_eq;
 use resolved::bodies::{Appearance, EARTH_EQUATORIAL_RADIUS, EARTH_POLAR_RADIUS, MOON_RADIUS};
 use resolved::gpu::bodies::{BodiesPass, BodiesUniforms, BodyRadiometry};
 use resolved::gpu::targets::FrameTargets;
-use resolved::gpu::{AdapterChoice, GpuContext};
+use resolved::gpu::GpuContext;
 use resolved::nalgebra::{UnitQuaternion, Vector3};
 use resolved::radiometry::{Radiometry, AU, SUN_RADIUS};
 use resolved::{CameraConfig, Pose, Scene};
@@ -29,10 +29,6 @@ impl Frame {
     }
 }
 
-fn contexts() -> [GpuContext; 2] {
-    [AdapterChoice::Gpu, AdapterChoice::Cpu].map(|c| GpuContext::new(&c).unwrap())
-}
-
 fn body_radiometry() -> BodyRadiometry {
     let radiometry = Radiometry::new(&CameraConfig::example());
     BodyRadiometry::new(&radiometry, &Appearance::default()).unwrap()
@@ -49,7 +45,12 @@ fn render_with(
     camera: &Pose,
     config: &CameraConfig,
 ) -> Frame {
-    let targets = FrameTargets::new(ctx, config.width, config.height);
+    let targets = FrameTargets::new(
+        ctx,
+        config.width,
+        config.height,
+        config.channels.len() as u32,
+    );
     let pass = BodiesPass::new(ctx);
     let uniforms = BodiesUniforms::new(
         scene,
@@ -100,8 +101,8 @@ fn sunlit_earth_center_radiance_and_depth() {
     let scene = scene(Vector3::new(AU, 0.0, 0.0), moon_away());
     let expected = body_radiometry().earth_lambert_1au;
     let near = CameraConfig::example().near_clip;
-    for ctx in contexts() {
+    for ctx in common::contexts() {
-        let frame = render(&ctx, "earth-sunlit", &scene, &earth_camera());
+        let frame = render(ctx, "earth-sunlit", &scene, &earth_camera());
         let (radiance, depth) = frame.at(960, 540);
         for k in 0..3 {
             assert_relative_eq!(radiance[k], expected[k], max_relative = 1e-4);
@@ -122,8 +123,8 @@ fn earth_limb_lands_where_projection_says() {
     let expected = (0..1920)
         .filter(|u| ((*u as f64 + 0.5) - k.cx).abs() < half_width)
         .count();
-    for ctx in contexts() {
+    for ctx in common::contexts() {
-        let frame = render(&ctx, "earth-limb", &scene, &earth_camera());
+        let frame = render(ctx, "earth-limb", &scene, &earth_camera());
         let hit = (0..frame.width)
             .filter(|&u| frame.at(u, 540).1 > 0.0)
             .count();
@@ -140,8 +141,8 @@ fn night_side_is_dark_and_hides_the_sun() {
     // The Sun directly behind Earth: the visible hemisphere is unlit and the solar disk
     // is eclipsed.
     let scene = scene(Vector3::new(-AU, 0.0, 0.0), moon_away());
-    for ctx in contexts() {
+    for ctx in common::contexts() {
-        let frame = render(&ctx, "earth-night", &scene, &earth_camera());
+        let frame = render(ctx, "earth-night", &scene, &earth_camera());
         assert!(
             frame.radiance.iter().all(|r| *r == [0.0; 4]),
             "{}",
@@ -162,8 +163,8 @@ fn solar_disk_size_and_radiance() {
     let radius = k.fx * (SUN_RADIUS / (sun - position).norm()).asin().tan();
     let expected = std::f64::consts::PI * radius * radius;
     let sun_radiance = body_radiometry().sun_radiance;
-    for ctx in contexts() {
+    for ctx in common::contexts() {
-        let frame = render(&ctx, "sun-disk", &scene, &camera);
+        let frame = render(ctx, "sun-disk", &scene, &camera);
         assert_eq!(frame.at(960, 540).0, sun_radiance);
         let lit = frame
             .radiance
@@ -186,8 +187,8 @@ fn moon_in_front_of_earth_wins() {
             .moon_lambert_1au
             .map(|l| l / (au * au) as f32)
     };
-    for ctx in contexts() {
+    for ctx in common::contexts() {
-        let frame = render(&ctx, "moon-in-front-of-earth", &scene, &earth_camera());
+        let frame = render(ctx, "moon-in-front-of-earth", &scene, &earth_camera());
         let (radiance, depth) = frame.at(960, 540);
         for k in 0..3 {
             assert_relative_eq!(radiance[k], expected[k], max_relative = 1e-4);
@@ -214,10 +215,10 @@ fn limb_follows_the_ellipsoid() {
     let x = a * a / RANGE;
     let north = Vector3::new(x, 0.0, b * (1.0 - (a / RANGE).powi(2)).sqrt());
     let east = Vector3::new(x, a * (1.0 - (a / RANGE).powi(2)).sqrt(), 0.0);
-    for ctx in contexts() {
+    for ctx in common::contexts() {
         // Boresight on the north limb, north up: Earth fills rows below the center.
         let camera = Pose::looking_at(position, north, Vector3::z()).unwrap();
-        let frame = render_with(&ctx, "limb-north", &scene, &camera, &config);
+        let frame = render_with(ctx, "limb-north", &scene, &camera, &config);
         let column: Vec<bool> = (0..1080).map(|v| frame.at(960, v).1 > 0.0).collect();
         let edge = column.iter().position(|hit| *hit).unwrap();
         assert!(column[edge..].iter().all(|hit| *hit));
@@ -229,7 +230,7 @@ fn limb_follows_the_ellipsoid() {
 
         // Boresight on the east limb, north up: Earth is to the left of center.
         let camera = Pose::looking_at(position, east, Vector3::z()).unwrap();
-        let frame = render_with(&ctx, "limb-east", &scene, &camera, &config);
+        let frame = render_with(ctx, "limb-east", &scene, &camera, &config);
         let row: Vec<bool> = (0..1920).map(|u| frame.at(u, 540).1 > 0.0).collect();
         let edge = row.iter().position(|hit| !*hit).unwrap();
         assert!(row[edge..].iter().all(|hit| !*hit));

crates/resolved/tests/common/mod.rs

@@ -1,14 +1,27 @@
-//! Shared by the GPU integration tests: saving what a test rendered, for looking at.
+//! Shared by the GPU integration tests: the devices to run on, and saving what a test
+//! rendered, for looking at.
 //!
 //! Set `RESOLVED_TEST_IMAGES=1` and each test writes PNGs to
 //! `target/tmp/test-images/`, one per test case and device.
 
-use std::path::PathBuf;
+// Each test binary compiles this module and uses only some of it.
+#![allow(dead_code)]
 
-use resolved::gpu::GpuContext;
+use std::path::PathBuf;
+use std::sync::OnceLock;
+
+use resolved::gpu::{AdapterChoice, GpuContext};
 use resolved::radiometry::{Radiometry, AU};
 use resolved::CameraConfig;
 
+/// A hardware GPU and lavapipe, opened once per test binary and shared by its tests.
+pub fn contexts() -> &'static [GpuContext; 2] {
+    static CONTEXTS: OnceLock<[GpuContext; 2]> = OnceLock::new();
+    CONTEXTS.get_or_init(|| {
+        [AdapterChoice::Gpu, AdapterChoice::Cpu].map(|c| GpuContext::new(&c).unwrap())
+    })
+}
+
 /// Writes `radiance` (example camera, row-major from the top) as an sRGB PNG named
 /// `<name>-<gpu|cpu>.png`, if `RESOLVED_TEST_IMAGES` is set.
 ///
@@ -19,32 +32,50 @@ pub fn save(name: &str, ctx: &GpuContext, radiance: &[[f32; 4]]) {
     if std::env::var_os("RESOLVED_TEST_IMAGES").is_none() {
         return;
     }
-    let config = CameraConfig::example();
+    let r = Radiometry::new(&CameraConfig::example());
-    let r = Radiometry::new(&config);
     let white: Vec<f32> = (0..3)
         .map(|k| r.lambert_radiance(k, 1.0, 1.0, AU) as f32)
         .collect();
-    let srgb = |x: f32| {
-        let x = x.clamp(0.0, 1.0);
-        let v = if x <= 0.003_130_8 {
-            12.92 * x
-        } else {
-            1.055 * x.powf(1.0 / 2.4) - 0.055
-        };
-        (v * 255.0).round() as u8
-    };
     let pixels: Vec<u8> = radiance
         .iter()
         .flat_map(|p| (0..3).map(move |k| (k, p[k])))
         .map(|(k, v)| srgb(v / white[k]))
         .collect();
+    write(name, ctx, &pixels);
+}
+
+/// Writes a DN image (example camera, three channels interleaved) as an sRGB PNG named
+/// `<name>-dn-<gpu|cpu>.png`, if `RESOLVED_TEST_IMAGES` is set. The ADC's full range
+/// maps to black through white, as a raw camera frame would display.
+pub fn save_dn(name: &str, ctx: &GpuContext, dn: &[u16]) {
+    if std::env::var_os("RESOLVED_TEST_IMAGES").is_none() {
+        return;
+    }
+    let max = CameraConfig::example().detector.dn_max() as f32;
+    let pixels: Vec<u8> = dn.iter().map(|&d| srgb(d as f32 / max)).collect();
+    write(&format!("{name}-dn"), ctx, &pixels);
+}
+
+/// Linear 0..1 to 8-bit sRGB.
+fn srgb(x: f32) -> u8 {
+    let x = x.clamp(0.0, 1.0);
+    let v = if x <= 0.003_130_8 {
+        12.92 * x
+    } else {
+        1.055 * x.powf(1.0 / 2.4) - 0.055
+    };
+    (v * 255.0).round() as u8
+}
+
+fn write(name: &str, ctx: &GpuContext, rgb: &[u8]) {
+    let config = CameraConfig::example();
     let dir = PathBuf::from(env!("CARGO_TARGET_TMPDIR")).join("test-images");
     std::fs::create_dir_all(&dir).unwrap();
     let device = if ctx.is_cpu() { "cpu" } else { "gpu" };
     let path = dir.join(format!("{name}-{device}.png"));
     image::save_buffer(
         &path,
-        &pixels,
+        rgb,
         config.width,
         config.height,
         image::ExtendedColorType::Rgb8,

crates/resolved/tests/mesh.rs

@@ -11,7 +11,7 @@ use resolved::eclipse::sun_visible_fraction;
 use resolved::gpu::bodies::{BodiesPass, BodiesUniforms, BodyRadiometry};
 use resolved::gpu::mesh::{MeshDraw, MeshPass};
 use resolved::gpu::targets::FrameTargets;
-use resolved::gpu::{AdapterChoice, GpuContext};
+use resolved::gpu::GpuContext;
 use resolved::model::{LoadOptions, Material, Mesh, Part};
 use resolved::nalgebra::{UnitQuaternion, Vector3};
 use resolved::radiometry::{Radiometry, AU};
@@ -30,10 +30,6 @@ impl Frame {
     }
 }
 
-fn contexts() -> [GpuContext; 2] {
-    [AdapterChoice::Gpu, AdapterChoice::Cpu].map(|c| GpuContext::new(&c).unwrap())
-}
-
 /// A cube with edges of 2 `scale` meters, centered on its origin, gray albedo 0.5.
 fn cube(scale: f64) -> Mesh {
     let mut s = String::from("solid cube\n");
@@ -88,7 +84,12 @@ fn render(ctx: &GpuContext, name: &str, scene: &Scene, camera: &Pose, mesh: &Mes
     let config = CameraConfig::example();
     let k = config.intrinsics();
     let radiometry = Radiometry::new(&config);
-    let targets = FrameTargets::new(ctx, config.width, config.height);
+    let targets = FrameTargets::new(
+        ctx,
+        config.width,
+        config.height,
+        config.channels.len() as u32,
+    );
     let bodies = BodiesPass::new(ctx);
     let mut meshes = MeshPass::new(ctx);
     let gpu_mesh = meshes.upload(ctx, mesh, &radiometry).unwrap();
@@ -142,8 +143,8 @@ fn sunlit_face_radiance_depth_and_silhouette() {
     let expected_width = (0..1920)
         .filter(|u| ((*u as f64 + 0.5) - k.cx).abs() < half_width)
         .count();
-    for ctx in contexts() {
+    for ctx in common::contexts() {
-        let frame = render(&ctx, "cube-sunlit", &scene, &camera, &cube(1.0));
+        let frame = render(ctx, "cube-sunlit", &scene, &camera, &cube(1.0));
         let (radiance, depth) = frame.at(960, 540);
         for c in 0..3 {
             assert_relative_eq!(radiance[c], expected[c], max_relative = 1e-4);
@@ -173,10 +174,10 @@ fn earth_hides_a_target_behind_it_but_not_in_front() {
     )
     .unwrap();
     let big = cube(5.0e4);
-    for ctx in contexts() {
+    for ctx in common::contexts() {
-        let empty = render(&ctx, "cube-scene-empty", &scene(sun, None), &camera, &big);
+        let empty = render(ctx, "cube-scene-empty", &scene(sun, None), &camera, &big);
         let behind = render(
-            &ctx,
+            ctx,
             "cube-behind-earth",
             &scene(sun, Some(identity_at(Vector3::new(-7.0e6, 0.0, 0.0)))),
             &camera,
@@ -188,7 +189,7 @@ fn earth_hides_a_target_behind_it_but_not_in_front() {
             ctx.info.name
         );
         let front = render(
-            &ctx,
+            ctx,
             "cube-in-front-of-earth",
             &scene(sun, Some(identity_at(Vector3::new(7.0e6, 0.0, 0.0)))),
             &camera,
@@ -210,14 +211,14 @@ fn penumbra_dims_and_umbra_darkens() {
     // Targets 7000 km behind Earth, one on the shadow axis and one at Earth's radius
     // off it (inside the penumbra), each seen from 30 m on its sunward side.
     let sun = Vector3::new(AU, 0.0, 0.0);
-    for ctx in contexts() {
+    for ctx in common::contexts() {
         for y in [0.0, EARTH_EQUATORIAL_RADIUS] {
             let target = Vector3::new(-7.0e6, y, 0.0);
             let camera =
                 Pose::looking_at(target + Vector3::new(30.0, 0.0, 0.0), target, Vector3::z())
                     .unwrap();
             let frame = render(
-                &ctx,
+                ctx,
                 &format!(
                     "cube-shadow-{}",
                     if y == 0.0 { "umbra" } else { "penumbra" }
@@ -271,9 +272,9 @@ fn thin_panel_is_lit_on_the_side_facing_the_camera() {
     let camera =
         Pose::looking_at(target + Vector3::new(30.0, 0.0, 0.0), target, Vector3::z()).unwrap();
     let expected = facing_sun((sun - target).norm());
-    for ctx in contexts() {
+    for ctx in common::contexts() {
         let frame = render(
-            &ctx,
+            ctx,
             "thin-panel",
             &scene(sun, Some(identity_at(target))),
             &camera,

crates/resolved/tests/sensor.rs

@@ -0,0 +1,166 @@
+//! The sensor pass against the CPU reference (`Radiometry::electrons` then
+//! `Detector::dn`), on the GPU and on the CPU (lavapipe).
+
+mod common;
+
+use resolved::bodies::Appearance;
+use resolved::gpu::bodies::{BodiesPass, BodiesUniforms, BodyRadiometry};
+use resolved::gpu::sensor::{SensorPass, SensorUniforms};
+use resolved::gpu::targets::FrameTargets;
+use resolved::gpu::GpuContext;
+use resolved::nalgebra::{UnitQuaternion, Vector2, Vector3};
+use resolved::radiometry::{Radiometry, AU};
+use resolved::{CameraConfig, Pose, Scene};
+
+/// Runs the sensor pass over `radiance` (example camera, rows from the top) and
+/// returns the DN image.
+fn sense(ctx: &GpuContext, name: &str, radiance: &[[f32; 4]]) -> Vec<u16> {
+    let config = CameraConfig::example();
+    let targets = FrameTargets::new(ctx, config.width, config.height, 3);
+    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(config.width * 16),
+            rows_per_image: Some(config.height),
+        },
+        targets.radiance.size(),
+    );
+    let pass = SensorPass::new(ctx);
+    let uniforms = SensorUniforms::new(&config, &Radiometry::new(&config));
+    let mut encoder = ctx.device.create_command_encoder(&Default::default());
+    pass.encode(ctx, &mut encoder, &targets, &uniforms);
+    ctx.queue.submit([encoder.finish()]);
+    let dn = targets.read_dn(ctx);
+    common::save_dn(name, ctx, &dn);
+    dn
+}
+
+/// The DN the CPU reference gives for `radiance` in channel `c` at pixel (u, v).
+fn reference(config: &CameraConfig, r: &Radiometry, radiance: f32, u: usize, v: usize) -> u16 {
+    let ray = config
+        .intrinsics()
+        .ray(&Vector2::new(u as f64 + 0.5, v as f64 + 0.5));
+    config.detector.dn(r.electrons(radiance as f64, ray.z))
+}
+
+/// Channel radiance of a white Lambertian surface facing the Sun at 1 AU.
+fn white() -> [f32; 4] {
+    let r = Radiometry::new(&CameraConfig::example());
+    let mut out = [0.0; 4];
+    for (k, o) in out.iter_mut().take(3).enumerate() {
+        *o = r.lambert_radiance(k, 1.0, 1.0, AU) as f32;
+    }
+    out
+}
+
+#[test]
+fn every_pixel_matches_the_cpu_reference() {
+    // A pattern from black to three times the white level (past saturation), different
+    // in each channel and pixel.
+    let config = CameraConfig::example();
+    let r = Radiometry::new(&config);
+    let w = white();
+    let level = |u: usize, v: usize, c: usize| ((u * 7 + v * 13 + c * 5) % 1000) as f32 / 333.0;
+    let radiance: Vec<[f32; 4]> = (0..1080)
+        .flat_map(|v| (0..1920).map(move |u| (u, v)))
+        .map(|(u, v)| [0, 1, 2, 3].map(|c| if c < 3 { w[c] * level(u, v, c) } else { 0.0 }))
+        .collect();
+    for ctx in common::contexts() {
+        let dn = sense(ctx, "sensor-pattern", &radiance);
+        let mut off_by_one = 0;
+        // Every 8th row, plus the last: enough to catch any indexing mistake.
+        for v in (0..1080).step_by(8).chain([1079]) {
+            for u in 0..1920 {
+                for c in 0..3 {
+                    let got = dn[(v * 1920 + u) * 3 + c];
+                    let want = reference(&config, &r, radiance[v * 1920 + u][c], u, v);
+                    assert!(
+                        got.abs_diff(want) <= 1,
+                        "{} ({u},{v},{c}): {got} vs {want}",
+                        ctx.info.name
+                    );
+                    off_by_one += usize::from(got != want);
+                }
+            }
+        }
+        // f32 on the GPU and f64 on the CPU disagree only where a value lands within
+        // rounding of a half-DN boundary: measured 557 of all 6.2 million values on the
+        // RTX 3060 and 699 on lavapipe, about 1 in 10,000.
+        assert!(
+            off_by_one < 1920 * 136 * 3 / 5000,
+            "{}: {off_by_one} off by one",
+            ctx.info.name
+        );
+    }
+}
+
+#[test]
+fn worked_example_white_surface_saturation_and_bias() {
+    let w = white();
+    for ctx in common::contexts() {
+        // A white surface facing the Sun fills the frame: docs/optics.md's worked example.
+        let dn = sense(ctx, "sensor-white", &vec![w; 1920 * 1080]);
+        let green = |u: usize, v: usize| dn[(v * 1920 + u) * 3 + 1];
+        assert_eq!(green(960, 540), 2292);
+        // cos^4 falloff: 6.9 degrees off axis at the corner, 10,819 electrons instead of
+        // 11,139, so 64 + 10,819 / 5 = 2227.8.
+        assert_eq!(green(0, 0), 2228);
+
+        // Looking straight at the Sun: the full well caps the signal at 20,000 / 5 + 64.
+        let sun = w.map(|x| x * 5e4);
+        assert!(sense(ctx, "sensor-sun", &vec![sun; 1920 * 1080])
+            .iter()
+            .all(|&d| d == 4064));
+        // Darkness reads as the bias.
+        assert!(sense(ctx, "sensor-dark", &vec![[0.0; 4]; 1920 * 1080])
+            .iter()
+            .all(|&d| d == 64));
+    }
+}
+
+#[test]
+fn sunlit_earth_through_the_whole_chain() {
+    // Earth seen from 200,000 km with the Sun behind the camera: the sub-camera point
+    // has albedo 0.3 facing the Sun, so 0.3 x 11,138.8 electrons in green: 732 DN.
+    let config = CameraConfig::example();
+    let radiometry = Radiometry::new(&config);
+    let scene = Scene {
+        earth_attitude: UnitQuaternion::identity(),
+        sun: Vector3::new(AU, 0.0, 0.0),
+        moon: Vector3::new(0.0, -3.844e8, 0.0),
+        targets: vec![],
+    };
+    let camera = Pose::looking_at(
+        Vector3::new(2.0e8, 0.0, 0.0),
+        Vector3::zeros(),
+        Vector3::z(),
+    )
+    .unwrap();
+    let body_radiometry = BodyRadiometry::new(&radiometry, &Appearance::default()).unwrap();
+    let bodies = BodiesUniforms::new(
+        &scene,
+        &camera,
+        &config.intrinsics(),
+        config.near_clip,
+        &body_radiometry,
+    );
+    let sensor = SensorUniforms::new(&config, &radiometry);
+    for ctx in common::contexts() {
+        let targets = FrameTargets::new(ctx, config.width, config.height, 3);
+        let mut encoder = ctx.device.create_command_encoder(&Default::default());
+        BodiesPass::new(ctx).encode(ctx, &mut encoder, &targets, &bodies);
+        SensorPass::new(ctx).encode(ctx, &mut encoder, &targets, &sensor);
+        ctx.queue.submit([encoder.finish()]);
+        let dn = targets.read_dn(ctx);
+        common::save_dn("sensor-earth", ctx, &dn);
+        assert_eq!(dn[(540 * 1920 + 960) * 3 + 1], 732);
+        assert_eq!(dn[0], 64, "space reads as the bias");
+    }
+}

docs/optics.md

@@ -242,7 +242,7 @@ the channel count generic for that reason.
 
 ## From radiance to electrons: the camera equation
 
-Code: `radiometry.rs` (`Radiometry::electrons`); the sensor pass, #9.
+Code: `radiometry.rs` (`Radiometry::electrons`); `gpu/sensor.rs`, `shaders/sensor.wgsl`.
 
 A lens with focal length f and aperture diameter D has f-number
 N = f / D. The example camera is f/4: a 50 mm lens with a 12.5 mm aperture.
@@ -282,7 +282,7 @@ 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`); the sensor pass, #9.
+Code: `camera.rs` (`Detector::dn`); `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.
@@ -294,7 +294,28 @@ Code: `camera.rs` (`Detector`); the sensor pass, #9.
   signal does not get clipped at 0. The example uses 64 DN.
 - **ADC bits.** The output resolution. 12 bits gives values 0 to 4095.
 
-So DN = bias + electrons / gain, rounded and clipped to the ADC range.
+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 (M3) will enter between the electrons and the gain.
+
+Two consequences of these numbers show up in every image:
+
+- **Saturation is not always the ADC's maximum.** In the example camera the
+  full well fills first: 20,000 / 5 + 64 = 4064 DN, below the 12-bit maximum
+  of 4095. A saturated pixel there reads 4064. Real cameras are designed
+  either way, and an algorithm that detects saturation by looking for 4095
+  would miss it in this one.
+- **Black is not zero.** A pixel that receives no light reads the bias,
+  64 DN, so empty space in a raw frame displays as dark gray.
+
+Raw frames are also not white-balanced. The camera reports what each
+channel collected, and the channels collect different amounts from the same
+gray surface. In the example camera, sunlight on a gray surface gives about
+15 % fewer electrons in blue than in red, so gray things look slightly
+yellow-green in a raw image. A white balance step, applied after capture,
+would scale the channels to make them equal. The renderer leaves it out,
+because it simulates what the sensor records.
 
 ### Worked example
 
@@ -308,6 +329,11 @@ checks the green figure.)
 - The full well is 20,000 electrons, so this surface is at 56 % of
   saturation.
 - Earth, with albedo near 0.3, reads about 64 + 3,342 / 5, about 732 DN.
+- The same white surface at the image corner, 6.9 degrees off axis, gets
+  cos^4 of that, 97 %: 10,819 electrons, 2228 DN.
+
+The sensor tests check the 2292, 2228, 732, 4064 and 64 DN figures
+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