PSF/Blurring. Turned out to be tough

crates/resolved-cli/src/main.rs

@@ -56,10 +56,14 @@ struct RenderArgs {
     /// needs Blue Marble data).
     #[arg(long, default_value = "L0")]
     lighting: LightingTier,
-    /// Camera fidelity tier: C0, or C1 (adds pixels that collect light over their area:
-    /// graded edges and continuous thin parts).
+    /// Camera fidelity tier: C0, or C1 (adds pixels that collect light over their area,
+    /// for graded edges and continuous thin parts, and lens blur).
     #[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<u32>,
     /// Device: auto, gpu, cpu, or part of an adapter's name.
     #[arg(long, default_value = "auto")]
     adapter: String,
@@ -193,6 +197,7 @@ fn render(args: &RenderArgs) -> Result<(), String> {
         profile: timings.is_some(),
         shadow_map_size: None,
         dn_in_place: None,
+        supersampling: args.supersampling,
     };
     let start = Instant::now();
     let mut renderer = Renderer::new(camera, options).map_err(|e| e.to_string())?;
@@ -254,11 +259,12 @@ fn render(args: &RenderArgs) -> Result<(), String> {
                 v[v.len() / 2]
             };
             println!(
-                "gpu      shadow {:.3} ms, earthshine {:.3} ms, bodies {:.3} ms, mesh {:.3} ms, sensor {:.3} ms, first to last {:.3} ms (medians)",
+                "gpu      shadow {:.3} ms, earthshine {:.3} ms, mesh {:.3} ms, bodies {:.3} ms, blur {:.3} ms, sensor {:.3} ms, first to last {:.3} ms (medians)",
                 median(|t| t.shadow),
                 median(|t| t.earthshine),
-                median(|t| t.bodies),
                 median(|t| t.mesh),
+                median(|t| t.bodies),
+                median(|t| t.blur),
                 median(|t| t.sensor),
                 median(|t| t.total)
             );

crates/resolved-cli/src/scene_file.rs

@@ -10,7 +10,8 @@ use std::path::{Path, PathBuf};
 use resolved::bodies::Appearance;
 use resolved::nalgebra::{Quaternion, UnitQuaternion, Vector2, Vector3};
 use resolved::{
-    CameraConfig, LoadOptions, Material, Mesh, ModelId, Pose, Renderer, Scene, Spectrum, Target,
+    CameraConfig, LoadOptions, Material, Mesh, ModelId, Pose, Psf, Renderer, Scene, Spectrum,
+    Target,
 };
 use serde::Deserialize;
 
@@ -60,6 +61,17 @@ pub struct CameraFile {
     pub bias: Option<f64>,
     pub full_well: Option<f64>,
     pub adc_bits: Option<u32>,
+    /// 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<Sigma>,
+}
+
+/// One sigma, or one per channel.
+#[derive(Debug, Deserialize)]
+#[serde(untagged)]
+pub enum Sigma {
+    All(f64),
+    PerChannel(Vec<f64>),
 }
 
 /// The observer's pose: `position`, plus either `attitude` or `look_at` (and `up`).
@@ -145,6 +157,16 @@ impl SceneFile {
         if let Some(t) = c.transmission {
             camera.transmission = Spectrum::constant(t);
         }
+        camera.psf = match &c.psf_sigma {
+            None => camera.psf,
+            Some(Sigma::All(s)) if *s == 0.0 => Psf::None,
+            Some(Sigma::All(s)) => Psf::Gaussian {
+                sigma: vec![*s; camera.channels.len()],
+            },
+            Some(Sigma::PerChannel(sigma)) => Psf::Gaussian {
+                sigma: sigma.clone(),
+            },
+        };
         camera
     }
 
@@ -368,6 +390,26 @@ mod tests {
         assert_eq!((camera.width, camera.height), (640, 1080));
         assert_eq!(camera.exposure, 2e-3);
         assert_eq!(camera.detector.adc_bits, 10);
+
+        for (line, psf) in [
+            (
+                "psf_sigma = 1.5",
+                Psf::Gaussian {
+                    sigma: vec![1.5; 3],
+                },
+            ),
+            (
+                "psf_sigma = [0.8, 0.7, 0.6]",
+                Psf::Gaussian {
+                    sigma: vec![0.8, 0.7, 0.6],
+                },
+            ),
+            ("psf_sigma = 0", Psf::None),
+        ] {
+            let file: SceneFile =
+                toml::from_str(&format!("{MINIMAL}\n[camera]\n{line}\n")).unwrap();
+            assert_eq!(file.camera().psf, psf, "{line}");
+        }
     }
 
     #[test]

crates/resolved/src/camera.rs

@@ -51,6 +51,62 @@ impl Detector {
     }
 }
 
+/// How the optics spread light from one point over the focal plane: the point spread
+/// function (PSF). Applied from camera tier C1 up.
+#[derive(Clone, Debug, Default, PartialEq)]
+pub enum Psf {
+    /// None: every point lands at one point, as through a pinhole.
+    #[default]
+    None,
+    /// A Gaussian spot. The sum of diffraction, aberrations, defocus, charge diffusion
+    /// and jitter is often close to one.
+    Gaussian {
+        /// Standard deviation in pixels, one per channel. Blur usually grows toward the
+        /// red.
+        sigma: Vec<f64>,
+    },
+}
+
+/// The largest Gaussian PSF sigma, in pixels. At 4 x 4 finer pixels it reaches 240 of
+/// them, which the blur's workgroup tile is sized for (`blur.wgsl`). Wider blur is
+/// cheaper by FFT.
+pub const MAX_PSF_SIGMA: f64 = 15.0;
+
+impl Psf {
+    /// How far the PSF reaches in pixels `1 / scale` across: 4 sigma of the widest
+    /// channel, rounded up, beyond which a Gaussian holds less than 1e-4 of its light.
+    /// 0 for none.
+    pub fn radius(&self, scale: u32) -> u32 {
+        match self {
+            Psf::None => 0,
+            Psf::Gaussian { sigma } => {
+                (4.0 * scale as f64 * sigma.iter().fold(0.0f64, |a, &b| a.max(b))).ceil() as u32
+            }
+        }
+    }
+
+    /// The weights of channel `channel`'s blur along one axis over pixels `1 / scale`
+    /// across, for offsets `-radius(scale)` to `radius(scale)` of them: the Gaussian at
+    /// each offset, scaled to sum to 1. A 2D blur is this along rows, then along
+    /// columns. `[1]` for none.
+    ///
+    /// A Gaussian sampled at whole pixels keeps its variance within 1 % of sigma^2
+    /// from sigma = 0.65 of those pixels up, but loses it below: 2.3 % at 0.6, 14 % at
+    /// 0.5 and 92 % at 0.3, where the spot is narrower than a pixel. Finer pixels move
+    /// that limit down.
+    pub fn weights(&self, channel: usize, scale: u32) -> Vec<f64> {
+        let Psf::Gaussian { sigma } = self else {
+            return vec![1.0];
+        };
+        let (r, s) = (self.radius(scale) as i64, sigma[channel] * scale as f64);
+        let w: Vec<f64> = (-r..=r)
+            .map(|k| (-((k * k) as f64) / (2.0 * s * s)).exp())
+            .collect();
+        let sum: f64 = w.iter().sum();
+        w.iter().map(|x| x / sum).collect()
+    }
+}
+
 /// A pinhole camera with radiometric optics and a detector.
 #[derive(Clone, Debug, PartialEq)]
 pub struct CameraConfig {
@@ -70,6 +126,8 @@ pub struct CameraConfig {
     pub exposure: f64,
     /// Fraction of light the optics pass, excluding per-channel filters.
     pub transmission: Spectrum,
+    /// How the optics blur the image.
+    pub psf: Psf,
     /// The detector's color channels.
     pub channels: Vec<Channel>,
     /// Electrons to DN.
@@ -86,6 +144,9 @@ impl CameraConfig {
     /// - 1920 x 1080 pixels, 5.5 um pitch, 50 mm focal length, so about 12.1 x 6.8
     ///   degrees field of view and 110 urad per pixel.
     /// - f/4, flat 90 % optics transmission.
+    /// - A Gaussian blur with sigma 0.5 px in every channel: a sharp camera, not
+    ///   derived from its optics. Diffraction alone at f/4 gives an Airy core of sigma
+    ///   0.17 px.
     /// - 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.
@@ -106,6 +167,9 @@ impl CameraConfig {
             f_number: 4.0,
             exposure: 1e-4,
             transmission: Spectrum::constant(0.9),
+            psf: Psf::Gaussian {
+                sigma: vec![0.5; 3],
+            },
             channels: vec![
                 channel("r", 620e-9),
                 channel("g", 540e-9),
@@ -143,6 +207,14 @@ impl CameraConfig {
         if self.channels.is_empty() {
             return fail("at least one channel is required");
         }
+        if let Psf::Gaussian { sigma } = &self.psf {
+            if sigma.len() != self.channels.len() {
+                return fail("the PSF needs one sigma per channel");
+            }
+            if !sigma.iter().all(|&s| positive(s) && s <= MAX_PSF_SIGMA) {
+                return fail("PSF sigma must be positive and at most 15 pixels");
+            }
+        }
         let d = &self.detector;
         if !positive(d.gain) || !positive(d.full_well) || !d.bias.is_finite() || d.bias < 0.0 {
             return fail("gain and full well must be positive and bias nonnegative");
@@ -248,6 +320,45 @@ mod tests {
         let mut camera = CameraConfig::example();
         camera.detector.adc_bits = 17;
         assert!(camera.validate().is_err());
+
+        let mut camera = CameraConfig::example();
+        camera.psf = Psf::Gaussian {
+            sigma: vec![0.7; 2],
+        };
+        assert!(camera.validate().is_err());
+
+        let mut camera = CameraConfig::example();
+        camera.psf = Psf::Gaussian {
+            sigma: vec![0.7, 0.0, 0.7],
+        };
+        assert!(camera.validate().is_err());
+    }
+
+    #[test]
+    fn psf_weights_are_a_normalized_gaussian() {
+        let psf = Psf::Gaussian {
+            sigma: vec![0.7, 1.0, 1.3],
+        };
+        assert_eq!(psf.radius(1), 6);
+        assert_eq!(psf.radius(2), 11);
+        for scale in [1, 2] {
+            let r = psf.radius(scale) as f64;
+            for (k, s) in [0.7, 1.0, 1.3].into_iter().enumerate() {
+                let w = psf.weights(k, scale);
+                assert_eq!(w.len() as f64, 2.0 * r + 1.0);
+                assert_relative_eq!(w.iter().sum::<f64>(), 1.0, epsilon = 1e-15);
+                // In pixels 1 / scale across.
+                let variance: f64 = w
+                    .iter()
+                    .enumerate()
+                    .map(|(i, w)| w * (i as f64 - r).powi(2))
+                    .sum();
+                let sigma = s * scale as f64;
+                assert_relative_eq!(variance, sigma * sigma, max_relative = 3e-3);
+            }
+        }
+        assert_eq!(Psf::None.radius(2), 0);
+        assert_eq!(Psf::None.weights(0, 2), [1.0]);
     }
 
     #[test]

crates/resolved/src/fidelity.rs

@@ -39,10 +39,12 @@ pub enum CameraTier {
     #[default]
     C0,
     /// C0, plus each pixel collecting light over its whole area, so edges are graded and
-    /// features thinner than a pixel stay continuous: target meshes are covered at 4
-    /// points per pixel, and 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.
+    /// features thinner than a pixel stay continuous, and the lens's blur (the camera's
+    /// [`Psf`](crate::Psf)). The scene is rendered finer than the pixels, 2 x 2 by
+    /// default, blurred, then averaged into pixels. 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.
     C1,
 }
 

crates/resolved/src/gpu/blur.rs

@@ -0,0 +1,198 @@
+//! The pass that blurs the rendered radiance by the camera's point spread function and
+//! averages it into the image's pixels.
+//!
+//! It runs after the bodies pass and before the sensor pass, so light spreads before
+//! the pixels collect it and saturate. The radiance is rendered `scale` times finer
+//! than the image along each axis, so the blur acts on the scene before the pixels
+//! average it, rather than on pixel averages, which loses where in a pixel light fell.
+//!
+//! A Gaussian is separable: blurring every row, then every column, with the 1D weights
+//! of [`Psf::weights`] gives the 2D Gaussian blur. Averaging `scale` x `scale` rendered
+//! pixels into one is separable too, so each pass does both at once: along its axis an
+//! image pixel is the sum of 2 r + `scale` rendered pixels, weighted by the Gaussian
+//! summed over the `scale` pixels it covers. The rendered radiance has a margin of r
+//! pixels on each side, so the image's edges get the light from just outside them.
+//!
+//! Each workgroup makes up to 256 pixels of one row or column and first loads the
+//! texels they share into workgroup memory (`blur.wgsl`). On the RTX 3060, without finer
+//! rendering, that was no faster than loading each tap from the texture at sigma
+//! 0.7 px, where both move the image through memory at the card's bandwidth, and 1.5
+//! to 1.7 times faster at 2 to 4 px.
+
+use super::profile::PassTimestamps;
+use super::targets::FrameTargets;
+use super::GpuContext;
+use crate::camera::Psf;
+
+/// Image pixels along the axis per workgroup at `scale` (`blur.wgsl`): as many as fit
+/// 256 rendered pixels, to share each loaded texel among as many as the tile allows.
+fn line(scale: u32) -> u32 {
+    match scale {
+        1 => 256,
+        2 => 128,
+        _ => 64,
+    }
+}
+
+/// The shader's uniform block; the layout matches `Params` in `blur.wgsl`.
+#[repr(C)]
+#[derive(Clone, Copy, Debug, PartialEq, bytemuck::Pod, bytemuck::Zeroable)]
+struct Params {
+    scale: u32,
+    radius: u32,
+    margin: u32,
+    unused: u32,
+}
+
+/// The compiled rows and columns pipelines and the weights.
+#[derive(Debug)]
+pub struct BlurPass {
+    line: u32,
+    rows: wgpu::ComputePipeline,
+    columns: wgpu::ComputePipeline,
+    weights: wgpu::Buffer,
+    params: wgpu::Buffer,
+}
+
+impl BlurPass {
+    /// Compiles the shader and uploads the weights of `psf` for `channels` channels,
+    /// over radiance rendered `scale` times finer than the image, with a margin of
+    /// `psf.radius(scale)`. `None` with no PSF and a `scale` of 1, where there is
+    /// nothing to do.
+    pub fn new(ctx: &GpuContext, psf: &Psf, channels: usize, scale: u32) -> Option<Self> {
+        use wgpu::util::DeviceExt;
+        if *psf == Psf::None && scale == 1 {
+            return None;
+        }
+        let per_channel: Vec<Vec<f64>> = (0..channels)
+            .map(|k| averaged(&psf.weights(k, scale), scale))
+            .collect();
+        let weights: Vec<[f32; 4]> = (0..per_channel[0].len())
+            .map(|i| std::array::from_fn(|k| per_channel.get(k).map_or(0.0, |w| w[i] as f32)))
+            .collect();
+        let radius = psf.radius(scale);
+        let params = Params {
+            scale,
+            radius,
+            margin: radius,
+            unused: 0,
+        };
+        let module = ctx
+            .device
+            .create_shader_module(wgpu::ShaderModuleDescriptor {
+                label: Some("blur"),
+                source: wgpu::ShaderSource::Wgsl(include_str!("../shaders/blur.wgsl").into()),
+            });
+        let pipeline = |entry| {
+            ctx.device
+                .create_compute_pipeline(&wgpu::ComputePipelineDescriptor {
+                    label: Some(entry),
+                    layout: None,
+                    module: &module,
+                    entry_point: Some(entry),
+                    compilation_options: wgpu::PipelineCompilationOptions {
+                        constants: &[("LINE", line(scale) as f64)],
+                        ..Default::default()
+                    },
+                    cache: None,
+                })
+        };
+        let buffer = |label, contents: &[u8], usage| {
+            ctx.device
+                .create_buffer_init(&wgpu::util::BufferInitDescriptor {
+                    label: Some(label),
+                    contents,
+                    usage,
+                })
+        };
+        Some(Self {
+            line: line(scale),
+            rows: pipeline("rows"),
+            columns: pipeline("columns"),
+            weights: buffer(
+                "blur weights",
+                bytemuck::cast_slice(&weights),
+                wgpu::BufferUsages::STORAGE,
+            ),
+            params: buffer(
+                "blur params",
+                bytemuck::bytes_of(&params),
+                wgpu::BufferUsages::UNIFORM,
+            ),
+        })
+    }
+
+    /// Records blurring and averaging `targets`' radiance into its blurred image. The
+    /// targets need the scale and margin this pass was made for.
+    pub fn encode(
+        &self,
+        ctx: &GpuContext,
+        encoder: &mut wgpu::CommandEncoder,
+        targets: &FrameTargets,
+        timestamps: Option<PassTimestamps>,
+    ) {
+        let [rows, image] = targets
+            .blurred
+            .as_ref()
+            .expect("targets allocated with a margin or a scale");
+        let mut pass = encoder.begin_compute_pass(&wgpu::ComputePassDescriptor {
+            label: Some("blur"),
+            timestamp_writes: timestamps.map(PassTimestamps::compute),
+        });
+        let (w, h, line) = (image.width(), rows.height(), self.line);
+        for (pipeline, source, destination, groups) in [
+            (&self.rows, &targets.radiance, rows, (w.div_ceil(line), h)),
+            (
+                &self.columns,
+                rows,
+                image,
+                (w, image.height().div_ceil(line)),
+            ),
+        ] {
+            let bind_group = ctx.device.create_bind_group(&wgpu::BindGroupDescriptor {
+                label: Some("blur"),
+                layout: &pipeline.get_bind_group_layout(0),
+                entries: &[
+                    wgpu::BindGroupEntry {
+                        binding: 0,
+                        resource: self.weights.as_entire_binding(),
+                    },
+                    wgpu::BindGroupEntry {
+                        binding: 1,
+                        resource: wgpu::BindingResource::TextureView(
+                            &source.create_view(&Default::default()),
+                        ),
+                    },
+                    wgpu::BindGroupEntry {
+                        binding: 2,
+                        resource: wgpu::BindingResource::TextureView(
+                            &destination.create_view(&Default::default()),
+                        ),
+                    },
+                    wgpu::BindGroupEntry {
+                        binding: 3,
+                        resource: self.params.as_entire_binding(),
+                    },
+                ],
+            });
+            pass.set_pipeline(pipeline);
+            pass.set_bind_group(0, &bind_group, &[]);
+            pass.dispatch_workgroups(groups.0, groups.1, 1);
+        }
+    }
+}
+
+/// The weights of one image pixel's taps along an axis: `weights` (2 r + 1 of them,
+/// over rendered pixels) summed over the `scale` rendered pixels the image pixel
+/// covers, and divided by `scale` to average them. 2 r + `scale` taps.
+pub fn averaged(weights: &[f64], scale: u32) -> Vec<f64> {
+    let s = scale as usize;
+    (0..weights.len() + s - 1)
+        .map(|t| {
+            (0..s)
+                .filter_map(|a| t.checked_sub(a).and_then(|i| weights.get(i)))
+                .sum::<f64>()
+                / s as f64
+        })
+        .collect()
+}

crates/resolved/src/gpu/bodies.rs

@@ -257,8 +257,25 @@ pub struct BodiesPass {
     bind_group: wgpu::BindGroup,
 }
 
+/// Where the samples of a pixel lie, from its top-left corner, in pixels: the standard
+/// patterns of Vulkan, Direct3D and Metal, which WebGPU requires for 4 samples. The
+/// mesh-edge tests check them on each device.
+fn sample_positions(samples: u32) -> &'static [[f32; 2]] {
+    match samples {
+        1 => &[[0.5, 0.5]],
+        2 => &[[0.75, 0.75], [0.25, 0.25]],
+        4 => &[
+            [0.375, 0.125],
+            [0.875, 0.375],
+            [0.125, 0.625],
+            [0.625, 0.875],
+        ],
+        _ => panic!("no sample pattern for {samples} samples"),
+    }
+}
+
 impl BodiesPass {
-    /// Compiles the shader for mesh targets of `samples` samples per pixel, casting grids
+    /// Compiles the shader for mesh targets of `samples` (1, 2 or 4) samples per pixel, casting grids
     /// of rays over the pixels on bodies' edges if `area`, and allocates the uniform
     /// buffer.
     pub fn new(ctx: &GpuContext, samples: u32, area: bool) -> Self {
@@ -271,9 +288,15 @@ impl BodiesPass {
         } else {
             ("texture_2d<f32>", "texture_depth_2d")
         };
+        let positions = sample_positions(samples)
+            .iter()
+            .map(|[x, y]| format!("vec2<f32>({x}, {y})"))
+            .collect::<Vec<_>>()
+            .join(", ");
         let prelude = format!(
             "const SAMPLES: u32 = {samples}u;\nconst AREA: bool = {area};\n\
-             alias MeshRadiance = {radiance};\nalias MeshDepth = {depth};\n"
+             alias MeshRadiance = {radiance};\nalias MeshDepth = {depth};\n\
+             const SAMPLE_POSITIONS = array<vec2<f32>, {samples}>({positions});\n"
         );
         let module = device.create_shader_module(wgpu::ShaderModuleDescriptor {
             label: Some("bodies"),

crates/resolved/src/gpu/mod.rs

@@ -1,5 +1,6 @@
 //! GPU plumbing: adapter selection, device setup, render targets and passes.
 
+pub mod blur;
 pub mod bodies;
 pub mod context;
 pub mod earthshine;

crates/resolved/src/gpu/profile.rs

@@ -8,7 +8,7 @@
 use super::GpuContext;
 
 /// The passes timed, in the order they run.
-pub const PASSES: [&str; 5] = ["shadow", "earthshine", "mesh", "bodies", "sensor"];
+pub const PASSES: [&str; 6] = ["shadow", "earthshine", "mesh", "bodies", "blur", "sensor"];
 
 /// Where one pass writes its begin and end timestamps.
 #[derive(Clone, Copy, Debug)]
@@ -71,6 +71,9 @@ pub struct PassTimings {
     pub bodies: f64,
     /// Target meshes.
     pub mesh: f64,
+    /// The point spread function and averaging finer rendering into pixels (zero work
+    /// below camera C1).
+    pub blur: f64,
     /// Radiance to DN.
     pub sensor: f64,
     /// From the start of the first pass to the end of the last, gaps included.
@@ -151,8 +154,9 @@ impl Profiler {
             earthshine: ms(ticks[2], ticks[3]),
             mesh: ms(ticks[4], ticks[5]),
             bodies: ms(ticks[6], ticks[7]),
-            sensor: ms(ticks[8], ticks[9]),
-            total: ms(ticks[0], ticks[9]),
+            blur: ms(ticks[8], ticks[9]),
+            sensor: ms(ticks[10], ticks[11]),
+            total: ms(ticks[0], ticks[11]),
         }
     }
 }

crates/resolved/src/gpu/sensor.rs

@@ -1,6 +1,6 @@
 //! 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
+//! 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`].
 
@@ -93,7 +93,7 @@ impl SensorPass {
     ) {
         ctx.queue
             .write_buffer(&self.uniforms, 0, bytemuck::bytes_of(uniforms));
-        let radiance = targets.radiance.create_view(&Default::default());
+        let radiance = targets.sensed().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),

crates/resolved/src/gpu/targets.rs

@@ -17,16 +17,27 @@ pub const DEPTH_COMPARE: wgpu::CompareFunction = wgpu::CompareFunction::Greater;
 /// Target meshes are drawn first, into their own radiance and depth with
 /// [`FrameTargets::samples`] samples per pixel. The bodies pass then draws Earth, the
 /// Moon and the Sun and combines them with the meshes' samples into one value per
-/// pixel, in `radiance` and `depth`, which the sensor pass reads.
+/// pixel, in `radiance` and `depth`. These are rendered [`FrameTargets::scale`] times
+/// finer than the image along each axis, and with a blur they extend
+/// [`FrameTargets::margin`] of those finer pixels past the image on every side, so that
+/// light from just outside the frame can blur into it. The blur passes then blur them
+/// and average them into the image's pixels, in `blurred`. The sensor pass reads
+/// [`FrameTargets::sensed`].
 #[derive(Debug)]
 pub struct FrameTargets {
     /// Image width in pixels.
     pub width: u32,
     /// Image height in pixels.
     pub height: u32,
-    /// Samples per pixel of the mesh targets: 1, or 4 to cover meshes at 4 points.
+    /// Rendered pixels along each axis per image pixel: 1, or more to render finer.
+    pub scale: u32,
+    /// Rendered pixels past the image on each side, for the blur; 0 without one.
+    pub margin: u32,
+    /// Samples per rendered pixel of the mesh targets: 1, or 2 or 4 to cover meshes at
+    /// several points.
     pub samples: u32,
-    /// The meshes' channel radiance, [`RADIANCE_FORMAT`], `samples` per pixel.
+    /// The meshes' channel radiance, [`RADIANCE_FORMAT`], `samples` per pixel. This and
+    /// the next three are rendered at `scale` and include the margin.
     pub mesh_radiance: wgpu::Texture,
     /// The meshes' depth, [`DEPTH_FORMAT`], `samples` per pixel.
     pub mesh_depth: wgpu::Texture,
@@ -35,6 +46,10 @@ pub struct FrameTargets {
     /// Depth of the nearest surface at each pixel, [`DEPTH_FORMAT`]: the body at the
     /// pixel's center or the nearest mesh sample in front of it.
     pub depth: wgpu::Texture,
+    /// With a blur or a `scale` above 1: `radiance` blurred and averaged along rows,
+    /// leaving the rows rendered, then along columns into the image,
+    /// [`RADIANCE_FORMAT`].
+    pub blurred: Option<[wgpu::Texture; 2]>,
     /// Number of camera channels.
     pub channels: u32,
     /// The sensor pass's output: u16 DN, channels interleaved per pixel, rows from the
@@ -43,10 +58,21 @@ pub struct FrameTargets {
 }
 
 impl FrameTargets {
-    /// Allocates targets of `width` x `height` pixels with `channels` channels, and
-    /// `samples` samples per pixel for the meshes.
-    pub fn new(ctx: &GpuContext, width: u32, height: u32, channels: u32, samples: u32) -> Self {
-        let texture = |label, format, usage, sample_count| {
+    /// Allocates targets for a `width` x `height` image with `channels` channels,
+    /// rendered `scale` times finer along each axis with `samples` samples per rendered
+    /// pixel for the meshes. A nonzero `margin` renders that many pixels past the image
+    /// on each side. The blur's images are allocated with a margin or a `scale` above 1.
+    pub fn new(
+        ctx: &GpuContext,
+        width: u32,
+        height: u32,
+        channels: u32,
+        samples: u32,
+        scale: u32,
+        margin: u32,
+    ) -> Self {
+        let (outer_width, outer_height) = (width * scale + 2 * margin, height * scale + 2 * margin);
+        let sized = |label, format, usage, sample_count, width, height| {
             ctx.device.create_texture(&wgpu::TextureDescriptor {
                 label: Some(label),
                 size: wgpu::Extent3d {
@@ -62,6 +88,9 @@ impl FrameTargets {
                 view_formats: &[],
             })
         };
+        let texture = |label, format, usage, samples| {
+            sized(label, format, usage, samples, outer_width, outer_height)
+        };
         let usage = wgpu::TextureUsages::RENDER_ATTACHMENT
             | wgpu::TextureUsages::TEXTURE_BINDING
             | wgpu::TextureUsages::COPY_SRC;
@@ -74,9 +103,23 @@ impl FrameTargets {
             } else {
                 wgpu::TextureUsages::empty()
             };
+        let blur = |label, height| {
+            sized(
+                label,
+                RADIANCE_FORMAT,
+                wgpu::TextureUsages::STORAGE_BINDING
+                    | wgpu::TextureUsages::TEXTURE_BINDING
+                    | wgpu::TextureUsages::COPY_SRC,
+                1,
+                width,
+                height,
+            )
+        };
         Self {
             width,
             height,
+            scale,
+            margin,
             samples,
             mesh_radiance: texture("mesh radiance", RADIANCE_FORMAT, mesh_usage, samples),
             mesh_depth: texture("mesh depth", DEPTH_FORMAT, mesh_usage, samples),
@@ -88,6 +131,8 @@ impl FrameTargets {
                 1,
             ),
             depth: texture("depth", DEPTH_FORMAT, usage, 1),
+            blurred: (margin > 0 || scale > 1)
+                .then(|| [blur("blurred rows", outer_height), blur("blurred", height)]),
             channels,
             dn: ctx.device.create_buffer(&wgpu::BufferDescriptor {
                 label: Some("dn"),
@@ -119,6 +164,15 @@ impl FrameTargets {
             .contains(wgpu::Features::MAPPABLE_PRIMARY_BUFFERS)
     }
 
+    /// The radiance image the sensor pass reads, `width` x `height`: the blurred image,
+    /// or without a blur, `radiance`.
+    pub fn sensed(&self) -> &wgpu::Texture {
+        match &self.blurred {
+            Some([_, image]) => image,
+            None => &self.radiance,
+        }
+    }
+
     /// 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
@@ -149,14 +203,14 @@ impl FrameTargets {
         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.
+    /// Reads the radiance image the sensor sees ([`FrameTargets::sensed`]) 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]> {
-        read_texture(ctx, &self.radiance, wgpu::TextureAspect::All, 16)
+        read_texture(ctx, self.sensed(), wgpu::TextureAspect::All, 16)
     }
 
-    /// Reads the depth target back, row by row from the top. Slow, like
-    /// [`FrameTargets::read_radiance`].
+    /// Reads the depth target back, as rendered (at `scale`, with the margin), row by
+    /// row from the top. Slow, like [`FrameTargets::read_radiance`].
     pub fn read_depth(&self, ctx: &GpuContext) -> Vec<f32> {
         read_texture(ctx, &self.depth, wgpu::TextureAspect::DepthOnly, 4)
     }

crates/resolved/src/lib.rs

@@ -37,7 +37,7 @@ mod solar_data;
 pub mod spectrum;
 pub mod sun;
 
-pub use camera::{CameraConfig, Channel, Detector, Intrinsics};
+pub use camera::{CameraConfig, Channel, Detector, Intrinsics, Psf};
 pub use earth::EarthTexture;
 pub use error::Error;
 pub use fidelity::{CameraTier, Fidelity, LightingTier};

crates/resolved/src/render.rs

@@ -12,12 +12,13 @@ use std::sync::mpsc;
 use nalgebra::Vector3;
 
 use crate::bodies::{Appearance, MOON_RADIUS};
-use crate::camera::CameraConfig;
+use crate::camera::{CameraConfig, Intrinsics, Psf};
 use std::sync::Arc;
 
 use crate::earth::EarthTexture;
 use crate::eclipse::{sun_visible_fraction, sun_visible_light};
 use crate::fidelity::{CameraTier, Fidelity, LightingTier};
+use crate::gpu::blur::BlurPass;
 use crate::gpu::bodies::{BodiesPass, BodiesUniforms, BodyRadiometry};
 use crate::gpu::earthshine::{EarthshinePass, TargetLight};
 use crate::gpu::mesh::{GpuMesh, MeshDraw, MeshPass, ShadowDraw};
@@ -47,8 +48,9 @@ pub struct RendererOptions {
     pub moon_texture: Option<Arc<MoonTexture>>,
     /// How Earth and the Moon reflect light.
     pub appearance: Appearance,
-    /// Also read back the channel radiance image each frame. Costs another 16 bytes per
-    /// pixel of transfer; for debugging and tests.
+    /// Also read back the channel radiance image each frame, as the pixels receive it
+    /// (blurred, from camera C1). Costs another 16 bytes per pixel of transfer; for
+    /// debugging and tests.
     pub read_radiance: bool,
     /// Time each pass on the GPU and report it in [`Frame::timings`]. Needs an adapter
     /// with timestamp queries; without one, timings are `None`.
@@ -64,6 +66,11 @@ pub struct RendererOptions {
     /// the RTX 3060, 1.9 ms on lavapipe), but wgpu warns it can be slower on GPUs with
     /// their own memory; this is the way to check, or to opt out.
     pub dn_in_place: Option<bool>,
+    /// From camera C1: render each pixel as n x n finer pixels, blur them by the camera's
+    /// PSF and average them into the pixel, so the blur acts on where light fell within
+    /// each pixel. `None` for the tier's choice: 2 at C1. 1 to 4; each step multiplies
+    /// rendering's time and memory by about n^2. Must be 1 (or `None`) at C0.
+    pub supersampling: Option<u32>,
 }
 
 /// One rendered image.
@@ -107,6 +114,7 @@ pub struct Renderer {
     targets: FrameTargets,
     bodies: BodiesPass,
     meshes: MeshPass,
+    blur: Option<BlurPass>,
     shadows: ShadowPass,
     earthshine: EarthshinePass,
     sensor: SensorPass,
@@ -129,9 +137,48 @@ impl Renderer {
         let radiometry = Radiometry::new(&camera);
         let mut body_radiometry = BodyRadiometry::new(&radiometry, &options.appearance)?;
         body_radiometry.bodies_eclipsed = options.fidelity.lighting >= LightingTier::L1;
-        // From camera C1, meshes are covered at 4 points per pixel.
+        // From camera C1, the scene is rendered `scale` times finer than the pixels, with
+        // meshes covered at several points per rendered pixel and a margin as wide as the
+        // PSF reaches, then blurred and averaged into pixels.
         let area = options.fidelity.camera >= CameraTier::C1;
-        let samples = if area { 4 } else { 1 };
+        let scale = match (options.fidelity.camera, options.supersampling) {
+            (CameraTier::C0, None | Some(1)) => 1,
+            (CameraTier::C0, Some(n)) => {
+                return Err(Error::Unsupported(format!(
+                    "camera C0 renders one point per pixel, so it cannot supersample {n} x {n}"
+                )))
+            }
+            (_, None) => 2,
+            (_, Some(n @ 1..=4)) => n,
+            (_, Some(n)) => {
+                return Err(Error::Unsupported(format!(
+                    "supersampling must be 1 to 4 pixels across, not {n}"
+                )))
+            }
+        };
+        // Mesh samples per rendered pixel: 4, or 2 when rendering finer, which at 2 x 2
+        // still puts samples at 4 distinct offsets across and down each image pixel. A
+        // device that cannot render 2 (lavapipe) gets 4.
+        let two = [RADIANCE_FORMAT, DEPTH_FORMAT]
+            .iter()
+            .all(|&format| ctx.multisamples(format, 2));
+        let samples = match (area, scale) {
+            (false, _) => 1,
+            (true, 1) => 4,
+            (true, _) if two => 2,
+            (true, _) => 4,
+        };
+        let psf = if area { &camera.psf } else { &Psf::None };
+        let blur = BlurPass::new(&ctx, psf, camera.channels.len(), scale);
+        let margin = psf.radius(scale);
+        let rendered = [camera.width, camera.height].map(|n| n * scale + 2 * margin);
+        let largest = ctx.device.limits().max_texture_dimension_2d;
+        if rendered.iter().any(|&n| n > largest) {
+            return Err(Error::Unsupported(format!(
+                "{} x {} pixels at {scale} x {scale} supersampling renders {} x {}, larger than {} can ({largest})",
+                camera.width, camera.height, rendered[0], rendered[1], ctx.info.name
+            )));
+        }
         for format in [RADIANCE_FORMAT, DEPTH_FORMAT] {
             if !ctx.multisamples(format, samples) {
                 return Err(Error::Unsupported(format!(
@@ -140,13 +187,24 @@ impl Renderer {
                 )));
             }
         }
+        // Rendering finer multiplies the targets' memory by scale^2: catch running out
+        // here rather than as wgpu's panic.
+        let scope = ctx.device.push_error_scope(wgpu::ErrorFilter::OutOfMemory);
         let targets = FrameTargets::new(
             &ctx,
             camera.width,
             camera.height,
             camera.channels.len() as u32,
             samples,
+            scale,
+            margin,
         );
+        if let Some(e) = pollster::block_on(scope.pop()) {
+            return Err(Error::Gpu(format!(
+                "{} x {} pixels rendered at {scale} x {scale} ({} x {}) need more memory than {} has: {e}",
+                camera.width, camera.height, rendered[0], rendered[1], ctx.info.name
+            )));
+        }
         let staging = |label, size| {
             ctx.device.create_buffer(&wgpu::BufferDescriptor {
                 label: Some(label),
@@ -198,6 +256,7 @@ impl Renderer {
         let mut renderer = Self {
             bodies: BodiesPass::new(&ctx, samples, area),
             meshes: MeshPass::new(&ctx, samples),
+            blur,
             shadows: ShadowPass::new(&ctx, options.shadow_map_size.unwrap_or(4096)),
             earthshine: EarthshinePass::new(&ctx),
             sensor: SensorPass::new(&ctx),
@@ -326,7 +385,20 @@ impl Renderer {
             )));
         }
         let camera = observer.then(&self.camera.mount);
-        let intrinsics = self.camera.intrinsics();
+        // The rendered images: finer by the scale, margin included.
+        let intrinsics = {
+            let k = self.camera.intrinsics();
+            let (s, m) = (self.targets.scale, self.targets.margin);
+            let f = s as f64;
+            Intrinsics {
+                fx: k.fx * f,
+                fy: k.fy * f,
+                cx: k.cx * f + m as f64,
+                cy: k.cy * f + m as f64,
+                width: k.width * s + 2 * m,
+                height: k.height * s + 2 * m,
+            }
+        };
         let near = self.camera.near_clip;
         let occulters = [
             (Vector3::zeros(), crate::bodies::EARTH_EQUATORIAL_RADIUS),
@@ -435,12 +507,20 @@ impl Renderer {
         );
         self.bodies
             .encode(ctx, &mut encoder, &self.targets, &bodies, stamps(3));
+        match &self.blur {
+            Some(blur) => blur.encode(ctx, &mut encoder, &self.targets, stamps(4)),
+            // Nothing to blur: an empty pass still writes the timestamps.
+            None => drop(encoder.begin_compute_pass(&wgpu::ComputePassDescriptor {
+                label: Some("blur (none)"),
+                timestamp_writes: stamps(4).map(|t| t.compute()),
+            })),
+        }
         self.sensor.encode(
             ctx,
             &mut encoder,
             &self.targets,
             &self.sensor_uniforms,
-            stamps(4),
+            stamps(5),
         );
         if let Some(profiler) = &self.profiler {
             profiler.resolve(&mut encoder);
@@ -457,7 +537,7 @@ impl Renderer {
         if let Some((buffer, row)) = &self.radiance_staging {
             encoder.copy_texture_to_buffer(
                 wgpu::TexelCopyTextureInfo {
-                    texture: &self.targets.radiance,
+                    texture: self.targets.sensed(),
                     mip_level: 0,
                     origin: wgpu::Origin3d::ZERO,
                     aspect: wgpu::TextureAspect::All,
@@ -470,7 +550,7 @@ impl Renderer {
                         rows_per_image: Some(self.camera.height),
                     },
                 },
-                self.targets.radiance.size(),
+                self.targets.sensed().size(),
             );
         }
         let index = ctx.queue.submit([encoder.finish()]);

crates/resolved/src/shaders/blur.wgsl

@@ -0,0 +1,67 @@
+// The point spread function and the pixels' area: the rendered radiance, `scale`
+// times finer than the image, blurred and averaged into the image's pixels, along rows
+// then along columns, each a compute pass.
+//
+// Along its axis, each image pixel is the weighted sum of 2 radius + scale rendered
+// pixels: the Gaussian, summed over the scale rendered pixels the image pixel covers
+// (blur.rs). The source extends `margin` rendered pixels past the image on both sides.
+// The rows pass leaves the columns as rendered, margin included, for the columns pass.
+
+struct Params {
+    // Rendered pixels per image pixel along each axis.
+    scale: u32,
+    // How far the blur reaches, in rendered pixels.
+    radius: u32,
+    // Rendered pixels past the image on each side, at least the radius.
+    margin: u32,
+    unused: u32,
+}
+
+// Per tap, the weight of each channel.
+@group(0) @binding(0) var<storage, read> weights: array<vec4<f32>>;
+@group(0) @binding(1) var source: texture_2d<f32>;
+@group(0) @binding(2) var destination: texture_storage_2d<rgba32float, write>;
+@group(0) @binding(3) var<uniform> p: Params;
+
+// Destination pixels along the axis per workgroup, set per scale (blur.rs), and the
+// source texels they read: LINE x scale + 2 radius, at most 256 + 2 x 240
+// (MAX_PSF_SIGMA in camera.rs).
+override LINE: u32 = 64u;
+var<workgroup> tile: array<vec4<f32>, 768>;
+
+// The destination pixels from `start` to LINE along `axis`, one per invocation `i`:
+// loads the source texels they read into the tile once, then sums from it.
+fn blur(start: vec2<u32>, i: u32, axis: vec2<u32>) {
+    let size = textureDimensions(destination);
+    let source_size = textureDimensions(source);
+    let taps = 2u * p.radius + p.scale;
+    // The first source texel the line reads. Across the axis the source and the
+    // destination line up.
+    let along = dot(start, axis);
+    let first = start + axis * (along * (p.scale - 1u) + p.margin - p.radius);
+    let length = LINE * p.scale + 2u * p.radius;
+    for (var j = i; j < length; j += LINE) {
+        let at = min(first + j * axis, source_size - 1u);
+        tile[j] = textureLoad(source, at, 0);
+    }
+    workgroupBarrier();
+    let at = start + i * axis;
+    if at.x >= size.x || at.y >= size.y {
+        return;
+    }
+    var sum = vec4<f32>(0.0);
+    for (var t = 0u; t < taps; t++) {
+        sum += weights[t] * tile[i * p.scale + t];
+    }
+    textureStore(destination, at, sum);
+}
+
+@compute @workgroup_size(LINE)
+fn rows(@builtin(workgroup_id) group: vec3<u32>, @builtin(local_invocation_index) i: u32) {
+    blur(vec2<u32>(group.x * LINE, group.y), i, vec2<u32>(1u, 0u));
+}
+
+@compute @workgroup_size(LINE)
+fn columns(@builtin(workgroup_id) group: vec3<u32>, @builtin(local_invocation_index) i: u32) {
+    blur(vec2<u32>(group.x, group.y * LINE), i, vec2<u32>(0u, 1u));
+}

crates/resolved/src/shaders/bodies.wgsl

@@ -8,8 +8,9 @@
 // of the meshes' SAMPLES samples in front of the bodies there replaces 1 / SAMPLES of
 // the bodies' light. Depth is written as near / z.
 //
-// The code that includes this file first declares SAMPLES, AREA, and the types of the
-// mesh targets, MeshRadiance and MeshDepth (multisampled when SAMPLES > 1).
+// The code that includes this file first declares SAMPLES, AREA, the types of the mesh
+// targets, MeshRadiance and MeshDepth (multisampled when SAMPLES > 1), and
+// SAMPLE_POSITIONS, where in a pixel each sample lies (bodies.rs).
 
 struct Uniforms {
     // fx, fy, cx, cy in pixels.
@@ -65,14 +66,6 @@ struct Uniforms {
 @group(1) @binding(0) var mesh_radiance: MeshRadiance;
 @group(1) @binding(1) var mesh_depth: MeshDepth;
 
-// Where the samples of a 4-sample pixel lie, from its top-left corner: the standard
-// pattern, which WebGPU requires for 4 samples.
-const SAMPLE_POSITIONS = array<vec2<f32>, 4>(
-    vec2<f32>(0.375, 0.125),
-    vec2<f32>(0.875, 0.375),
-    vec2<f32>(0.125, 0.625),
-    vec2<f32>(0.625, 0.875),
-);
 // Nothing hit.
 const FAR: f32 = 3.4e38;
 
@@ -258,11 +251,18 @@ fn pixel_angle() -> f32 {
 }
 
 // Which of the SAMPLES mesh samples lies nearest `offset` from the pixel's top-left
-// corner, the first on a tie.
+// corner. On a tie, the first, except with 2 samples, whose ties lie along the
+// diagonal between them: there the left half of the diagonal goes to the upper-left
+// sample, so that each of the two is nearest to the same number of rays of a grid.
 fn nearest_sample(offset: vec2<f32>) -> u32 {
     if SAMPLES == 1u {
         return 0u;
     }
+    if SAMPLES == 2u {
+        let side = offset.x + offset.y - 1.0;
+        let upper_left = side < 0.0 || (side == 0.0 && offset.x < 0.5);
+        return select(0u, 1u, upper_left);
+    }
     var best = 0u;
     var best_distance = 2.0;
     for (var s = 0u; s < SAMPLES; s++) {
@@ -290,7 +290,7 @@ fn limb(c: vec3<f32>, r: f32, d: vec3<f32>) -> vec2<f32> {
 // sample it coarsely, so that at least 8 rays span its radius (to 64 at most). A grid
 // of 4 measured 45 % off in the total light of a disk 0.225 pixels in radius, 11 % at
 // 0.5 and 2.6 % at 1; 8 rays per radius keeps it under 1 % from 0.25 pixels up. Always
-// even, so that each of the 4 mesh samples is nearest to the same number of rays.
+// even, so that each mesh sample is nearest to the same number of rays.
 fn grid_size(d: vec3<f32>, pixel: f32) -> u32 {
     let earth = earth_stretched(d);
     var edges = array<vec2<f32>, 3>(

crates/resolved/tests/gpu/bodies.rs

@@ -64,6 +64,8 @@ fn render_limb(
         config.height,
         config.channels.len() as u32,
         1,
+        1,
+        0,
     );
     let pass = BodiesPass::new(ctx, 1, false);
     pass.set_limb_darkening(ctx, limb);

crates/resolved/tests/gpu/main.rs

@@ -13,6 +13,7 @@ mod gpu_context;
 mod mesh;
 mod moon_texture;
 mod pixel_area;
+mod psf;
 mod renderer;
 mod sensor;
 mod shadows;

crates/resolved/tests/gpu/mesh.rs

@@ -97,6 +97,8 @@ fn render(ctx: &GpuContext, name: &str, scene: &Scene, camera: &Pose, mesh: &Mes
         config.height,
         config.channels.len() as u32,
         1,
+        1,
+        0,
     );
     let bodies = BodiesPass::new(ctx, 1, false);
     let mut meshes = MeshPass::new(ctx, 1);

crates/resolved/tests/gpu/pixel_area.rs

@@ -11,8 +11,8 @@ use resolved::gpu::AdapterChoice;
 use resolved::nalgebra::{UnitQuaternion, Vector2, Vector3};
 use resolved::radiometry::{Radiometry, AU};
 use resolved::{
-    CameraConfig, CameraTier, Fidelity, Intrinsics, LoadOptions, Material, Mesh, Pose, Renderer,
-    RendererOptions, Scene, Target,
+    CameraConfig, CameraTier, Fidelity, Intrinsics, LoadOptions, Material, Mesh, Pose, Psf,
+    Renderer, RendererOptions, Scene, Target,
 };
 
 const ADAPTERS: [AdapterChoice; 2] = [AdapterChoice::Gpu, AdapterChoice::Cpu];
@@ -20,7 +20,20 @@ const ADAPTERS: [AdapterChoice; 2] = [AdapterChoice::Gpu, AdapterChoice::Cpu];
 /// The channel the tests compare (green).
 const CHANNEL: usize = 1;
 
+/// A renderer at camera `tier` with no blur, which would spread the edges these tests
+/// measure, and one rendered pixel per pixel (`psf.rs` tests both).
 fn renderer(adapter: &AdapterChoice, camera: CameraConfig, tier: CameraTier) -> Renderer {
+    renderer_at(adapter, camera, tier, 1)
+}
+
+/// [`renderer`], rendering `supersampling` x `supersampling` pixels per pixel.
+fn renderer_at(
+    adapter: &AdapterChoice,
+    mut camera: CameraConfig,
+    tier: CameraTier,
+    supersampling: u32,
+) -> Renderer {
+    camera.psf = Psf::None;
     let options = RendererOptions {
         adapter: adapter.clone(),
         fidelity: Fidelity {
@@ -28,6 +41,7 @@ fn renderer(adapter: &AdapterChoice, camera: CameraConfig, tier: CameraTier) ->
             ..Default::default()
         },
         read_radiance: true,
+        supersampling: Some(supersampling),
         ..Default::default()
     };
     Renderer::new(camera, options).unwrap()
@@ -83,40 +97,53 @@ fn scene(sun: Vector3<f64>, moon: Vector3<f64>, targets: Vec<Target>) -> Scene {
 #[test]
 fn a_mesh_edge_grades_with_the_covered_fraction() {
     // A plate's left edge swept across pixel column 960 in tenths of a pixel. With 4
-    // samples per pixel the pixel reads the covered fraction to within half of 1/4.
+    // samples per pixel, or 2 x 2 rendered pixels of 2 samples each, the samples lie at
+    // 4 distinct offsets across the pixel, and it reads the covered fraction to within
+    // half of 1/4. Lavapipe renders 4 samples at 2 x 2, at 8 offsets.
     let _devices = common::device_lock();
     let camera = deep_space();
     let k = CameraConfig::example().intrinsics();
     let range = 100.0;
     for adapter in &ADAPTERS {
-        let mut r = renderer(adapter, CameraConfig::example(), CameraTier::C1);
-        let model = r.load_model(&plate(1.0)).unwrap();
-        let name = r.adapter_info().name.clone();
-        for step in 0..10 {
-            let f = step as f64 / 10.0;
-            let edge = Vector3::new((960.0 + f - k.cx) / k.fx * range, 0.0, range);
-            let target = Target {
-                model,
-                pose: Pose::new(camera.position + edge, UnitQuaternion::identity()),
-            };
-            // The Sun behind the camera, lighting the plate face on.
-            let sun = camera.position - Vector3::z() * AU;
-            let moon = Vector3::new(0.0, -3.844e8, 0.0);
-            let image = render(&mut r, &scene(sun, moon, vec![target]), &camera);
-            let inside = image[540 * 1920 + 1000];
-            let edge_pixel = image[540 * 1920 + 960];
-            assert!(inside > 0.0, "{name}");
-            let fraction = edge_pixel / inside;
-            assert!(
-                (fraction - (1.0 - f)).abs() <= 0.125 + 1e-4,
-                "{name}: edge {f} px into the pixel: {fraction}, covered {}",
-                1.0 - f
+        for supersampling in [1, 2] {
+            let mut r = renderer_at(
+                adapter,
+                CameraConfig::example(),
+                CameraTier::C1,
+                supersampling,
             );
-            // A whole number of samples.
-            assert!(
-                (fraction * 4.0 - (fraction * 4.0).round()).abs() < 1e-4,
-                "{name}: {fraction}"
+            let model = r.load_model(&plate(1.0)).unwrap();
+            let name = format!(
+                "{}, {supersampling} x {supersampling}",
+                r.adapter_info().name
             );
+            for step in 0..10 {
+                let f = step as f64 / 10.0;
+                let edge = Vector3::new((960.0 + f - k.cx) / k.fx * range, 0.0, range);
+                let target = Target {
+                    model,
+                    pose: Pose::new(camera.position + edge, UnitQuaternion::identity()),
+                };
+                // The Sun behind the camera, lighting the plate face on.
+                let sun = camera.position - Vector3::z() * AU;
+                let moon = Vector3::new(0.0, -3.844e8, 0.0);
+                let image = render(&mut r, &scene(sun, moon, vec![target]), &camera);
+                let inside = image[540 * 1920 + 1000];
+                let edge_pixel = image[540 * 1920 + 960];
+                assert!(inside > 0.0, "{name}");
+                let fraction = edge_pixel / inside;
+                assert!(
+                    (fraction - (1.0 - f)).abs() <= 0.125 + 1e-4,
+                    "{name}: edge {f} px into the pixel: {fraction}, covered {}",
+                    1.0 - f
+                );
+                // A whole number of samples, of 4 or 2 x 2 x 4.
+                let n = (4 * supersampling * supersampling) as f64;
+                assert!(
+                    (fraction * n - (fraction * n).round()).abs() < 1e-4,
+                    "{name}: {fraction}"
+                );
+            }
         }
     }
 }

crates/resolved/tests/gpu/psf.rs

@@ -0,0 +1,398 @@
+//! Camera tier C1's blur by the point spread function, over radiance rendered finer
+//! than the pixels: the blur passes against the exact blurred and pixel-averaged scene
+//! and against a direct computation, and the renderer applying them, on the GPU and on
+//! the CPU (lavapipe).
+
+use crate::common;
+
+use resolved::camera::MAX_PSF_SIGMA;
+use resolved::gpu::blur::{averaged, BlurPass};
+use resolved::gpu::targets::FrameTargets;
+use resolved::gpu::{AdapterChoice, GpuContext};
+use resolved::nalgebra::{UnitQuaternion, Vector2, Vector3};
+use resolved::radiometry::AU;
+use resolved::{
+    CameraConfig, CameraTier, Error, Fidelity, Frame, Pose, Psf, Renderer, RendererOptions, Scene,
+};
+
+/// The image's size in the pass tests.
+const WIDTH: u32 = 48;
+const HEIGHT: u32 = 40;
+
+/// A sharp camera's blur, narrower than a pixel, the same in every channel.
+fn sharp() -> Psf {
+    Psf::Gaussian {
+        sigma: vec![0.4; 3],
+    }
+}
+
+/// Blurs `radiance`, rendered `scale` times finer than an image of `size` pixels with a
+/// margin of `psf.radius(scale)` (rows from the top), into the image.
+fn blur(
+    ctx: &GpuContext,
+    psf: &Psf,
+    scale: u32,
+    size: (u32, u32),
+    radiance: &[[f32; 4]],
+) -> Vec<[f32; 4]> {
+    let targets = FrameTargets::new(ctx, size.0, size.1, 3, 1, scale, psf.radius(scale));
+    let (w, h) = (targets.radiance.width(), targets.radiance.height());
+    assert_eq!(radiance.len(), (w * h) 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(w * 16),
+            rows_per_image: Some(h),
+        },
+        targets.radiance.size(),
+    );
+    let pass = BlurPass::new(ctx, psf, 3, scale).unwrap();
+    let mut encoder = ctx.device.create_command_encoder(&Default::default());
+    pass.encode(ctx, &mut encoder, &targets, None);
+    ctx.queue.submit([encoder.finish()]);
+    targets.read_radiance(ctx)
+}
+
+/// Rendered radiance `scale` times finer than the `WIDTH` x `HEIGHT` image, with
+/// `psf`'s margin: `f(x, y, d)` at each rendered pixel, given the image coordinates of
+/// its left and top edges and its width `d` in image pixels.
+fn render_fine(psf: &Psf, scale: u32, mut f: impl FnMut(f64, f64, f64) -> f64) -> Vec<[f32; 4]> {
+    let m = psf.radius(scale) as i64;
+    let (w, h) = (
+        (WIDTH * scale) as i64 + 2 * m,
+        (HEIGHT * scale) as i64 + 2 * m,
+    );
+    let d = 1.0 / scale as f64;
+    (0..w * h)
+        .map(|i| {
+            let (x, y) = (((i % w) - m) as f64 * d, ((i / w) - m) as f64 * d);
+            [f(x, y, d) as f32; 4]
+        })
+        .collect()
+}
+
+/// The standard normal cumulative distribution.
+fn phi(x: f64) -> f64 {
+    0.5 * (1.0 + erf(x / std::f64::consts::SQRT_2))
+}
+
+/// The error function, to 1.5e-7 (Abramowitz and Stegun 7.1.26).
+fn erf(x: f64) -> f64 {
+    let t = 1.0 / (1.0 + 0.3275911 * x.abs());
+    let poly = t
+        * (0.254829592
+            + t * (-0.284496736 + t * (1.421413741 + t * (-1.453152027 + t * 1.061405429))));
+    (1.0 - poly * (-x * x).exp()).copysign(x)
+}
+
+/// The fraction of a Gaussian of `sigma` centered at `x0` that lands in pixel
+/// `[u, u + 1)`.
+fn in_pixel(u: f64, x0: f64, sigma: f64) -> f64 {
+    phi((u + 1.0 - x0) / sigma) - phi((u - x0) / sigma)
+}
+
+#[test]
+fn points_land_as_the_blurred_spot_over_each_pixel() {
+    // Sixteen points of light 1 at positions across a pixel, each drawn into the
+    // rendered pixel it falls in. Each image pixel should get the share of a Gaussian
+    // of sigma 0.4 px centered on the point that lands in it. Rendering finer places
+    // each point more closely: the worst pixel's error, as a fraction of a point's
+    // light, measured 0.507, 0.191 and 0.045 at 1, 2 and 4 rendered pixels across.
+    let psf = sharp();
+    let points: Vec<(f64, f64)> = (0..16)
+        .map(|i| {
+            let (a, b) = ((i % 4) as f64, (i / 4) as f64);
+            (
+                6.0 + 11.0 * a + 0.13 + 0.25 * a,
+                6.0 + 9.0 * b + 0.07 + 0.25 * b,
+            )
+        })
+        .collect();
+    for ctx in common::contexts() {
+        let name = &ctx.info.name;
+        for (scale, tolerance) in [(1, 0.52), (2, 0.2), (4, 0.05)] {
+            let radiance = render_fine(&psf, scale, |x, y, d| {
+                let hit = points
+                    .iter()
+                    .any(|&(px, py)| (x..x + d).contains(&px) && (y..y + d).contains(&py));
+                if hit {
+                    1.0 / (d * d)
+                } else {
+                    0.0
+                }
+            });
+            let image = blur(ctx, &psf, scale, (WIDTH, HEIGHT), &radiance);
+            let mut worst = 0.0f64;
+            for (i, p) in image.iter().enumerate() {
+                let (u, v) = ((i % WIDTH as usize) as f64, (i / WIDTH as usize) as f64);
+                let expected: f64 = points
+                    .iter()
+                    .map(|&(px, py)| in_pixel(u, px, 0.4) * in_pixel(v, py, 0.4))
+                    .sum();
+                worst = worst.max((p[1] as f64 - expected).abs());
+            }
+            let total: f64 = image.iter().map(|p| p[1] as f64).sum();
+            assert!((total / 16.0 - 1.0).abs() < 1e-5, "{name}: {total}");
+            assert!(worst < tolerance, "{name}: {scale} x {scale}: {worst}");
+        }
+    }
+}
+
+#[test]
+fn an_edge_lands_as_the_blurred_edge_over_each_pixel() {
+    // A vertical edge, dark left of x0 and 1 to the right, with each rendered pixel
+    // holding the fraction of it right of x0, as C1 renders it. Each image pixel should
+    // be the edge blurred by a Gaussian of sigma 0.4 px and averaged over the pixel,
+    // sigma (G((u + 1 - x0) / sigma) - G((u - x0) / sigma)) with G(z) = z Phi(z) + phi(z).
+    // Measured in f64 at 129 edge positions: worst 0.118, 0.022 and 0.005 of the step
+    // at 1, 2 and 4 rendered pixels across, each with the edge on a pixel boundary.
+    let psf = sharp();
+    let sigma = 0.4;
+    let integral = |z: f64| z * phi(z) + (-z * z / 2.0).exp() / (2.0 * std::f64::consts::PI).sqrt();
+    for ctx in common::contexts() {
+        let name = &ctx.info.name;
+        for (scale, tolerance) in [(1, 0.119), (2, 0.0224), (4, 0.0052)] {
+            for x0 in [20.0, 20.1, 20.25, 20.3, 20.5, 20.875] {
+                let radiance =
+                    render_fine(&psf, scale, |x, _, d| ((x + d - x0) / d).clamp(0.0, 1.0));
+                let image = blur(ctx, &psf, scale, (WIDTH, HEIGHT), &radiance);
+                for v in [0, HEIGHT as usize / 2, HEIGHT as usize - 1] {
+                    for u in 10..30 {
+                        let x = u as f64;
+                        let expected =
+                            sigma * (integral((x + 1.0 - x0) / sigma) - integral((x - x0) / sigma));
+                        let got = image[v * WIDTH as usize + u][0] as f64;
+                        assert!(
+                            (got - expected).abs() < tolerance,
+                            "{name}: {scale} x {scale}, edge at {x0}, pixel ({u}, {v}): {got} vs {expected}"
+                        );
+                    }
+                }
+            }
+        }
+    }
+}
+
+#[test]
+fn separable_passes_match_a_direct_computation() {
+    // Rows then columns against blurring the rendered pixels by the 2D kernel
+    // w(i) w(j) and averaging each block into an image pixel, in f64, on random
+    // radiance: equal to 1e-6 of the brightest pixel, with a different sigma per
+    // channel and at 1, 2 and 3 rendered pixels across.
+    let psf = Psf::Gaussian {
+        sigma: vec![0.4, 0.7, 1.1],
+    };
+    let mut state = 0x2545_f491_4f6c_dd1du64;
+    let mut random = move || {
+        state ^= state << 13;
+        state ^= state >> 7;
+        state ^= state << 17;
+        (state >> 11) as f64 / (1u64 << 53) as f64
+    };
+    for scale in [1, 2, 3] {
+        let s = scale as usize;
+        let w = WIDTH as usize * s + 2 * psf.radius(scale) as usize;
+        let radiance = render_fine(&psf, scale, |_, _, _| random());
+        let peak = radiance.iter().flatten().fold(0.0f32, |a, &b| a.max(b)) as f64;
+        let weights: Vec<Vec<f64>> = (0..3).map(|k| psf.weights(k, scale)).collect();
+        let mut expected = vec![[0.0f64; 3]; (WIDTH * HEIGHT) as usize];
+        for (i, e) in expected.iter_mut().enumerate() {
+            let (u, v) = (i % WIDTH as usize, i / WIDTH as usize);
+            for (k, wk) in weights.iter().enumerate() {
+                let mut sum = 0.0;
+                // Rendered pixel (u s + a, v s + b) of the block, blurred.
+                for (a, b) in (0..s).flat_map(|a| (0..s).map(move |b| (a, b))) {
+                    for (j, wj) in wk.iter().enumerate() {
+                        for (i, wi) in wk.iter().enumerate() {
+                            let (x, y) = (u * s + a + i, v * s + b + j);
+                            sum += wj * wi * radiance[y * w + x][k] as f64;
+                        }
+                    }
+                }
+                e[k] = sum / (s * s) as f64;
+            }
+        }
+        for ctx in common::contexts() {
+            let name = &ctx.info.name;
+            let image = blur(ctx, &psf, scale, (WIDTH, HEIGHT), &radiance);
+            let worst = image
+                .iter()
+                .zip(&expected)
+                .flat_map(|(p, e)| (0..3).map(move |k| (p[k] as f64 - e[k]).abs()))
+                .fold(0.0f64, f64::max);
+            assert!(worst < 1e-6 * peak, "{name}: {scale} x {scale}: {worst}");
+        }
+    }
+}
+
+#[test]
+fn averaged_weights_cover_the_pixel() {
+    // 1, 2, 1 over half pixels: each image pixel takes the blur of two rendered pixels.
+    let w = averaged(&[0.25, 0.5, 0.25], 2);
+    assert_eq!(w, [0.125, 0.375, 0.375, 0.125]);
+    assert_eq!(averaged(&[1.0], 3), [1.0 / 3.0; 3]);
+}
+
+/// Blurs radiance rising linearly across an image of `size` pixels and its margin,
+/// rendered `scale` times finer, and checks the image has the plane's values at its
+/// pixels' centers: a symmetric blur that sums to 1, and the average over a pixel,
+/// leave a plane unchanged.
+fn check_plane(ctx: &GpuContext, psf: &Psf, scale: u32, size: (u32, u32)) {
+    let (s, m) = (scale as f64, psf.radius(scale) as f64);
+    let (w, h) = (
+        (size.0 * scale) as usize + 2 * m as usize,
+        (size.1 * scale) as usize + 2 * m as usize,
+    );
+    // In image pixels.
+    let plane = |x: f64, y: f64| 1.0 + 0.01 * x + 0.02 * y;
+    let radiance: Vec<[f32; 4]> = (0..w * h)
+        .map(|i| {
+            let x = ((i % w) as f64 - m + 0.5) / s;
+            let y = ((i / w) as f64 - m + 0.5) / s;
+            [plane(x, y) as f32; 4]
+        })
+        .collect();
+    let image = blur(ctx, psf, scale, size, &radiance);
+    for (i, p) in image.iter().enumerate() {
+        let (u, v) = (i % size.0 as usize, i / size.0 as usize);
+        let expected = plane(u as f64 + 0.5, v as f64 + 0.5);
+        for (k, &l) in p.iter().take(3).enumerate() {
+            assert!(
+                (l as f64 / expected - 1.0).abs() < 2e-6,
+                "{}: {psf:?} at {scale}, pixel ({u}, {v}), channel {k}: {l} vs {expected}",
+                ctx.info.name
+            );
+        }
+    }
+}
+
+#[test]
+fn light_from_the_margin_blurs_into_the_image() {
+    // A plane keeps its values out to the image's edges, which take light from the
+    // margin. Also with the widest blur at 4 x 4, whose taps fill the shader's tile on
+    // an image more than 64 pixels across.
+    let widest = Psf::Gaussian {
+        sigma: vec![MAX_PSF_SIGMA, 0.4, 1.0],
+    };
+    for ctx in common::contexts() {
+        for scale in [1, 2, 4] {
+            check_plane(ctx, &sharp(), scale, (WIDTH, HEIGHT));
+        }
+        check_plane(ctx, &widest, 4, (100, 70));
+    }
+}
+
+const ADAPTERS: [AdapterChoice; 2] = [AdapterChoice::Gpu, AdapterChoice::Cpu];
+
+fn renderer(
+    adapter: &AdapterChoice,
+    camera: CameraConfig,
+    tier: CameraTier,
+    supersampling: Option<u32>,
+) -> Result<Renderer, Error> {
+    let options = RendererOptions {
+        adapter: adapter.clone(),
+        fidelity: Fidelity {
+            camera: tier,
+            ..Default::default()
+        },
+        read_radiance: true,
+        supersampling,
+        ..Default::default()
+    };
+    Renderer::new(camera, options)
+}
+
+/// A full Moon 41 pixels in radius, centered 6 pixels right of and 4 below the image
+/// center, seen from 4.2e4 km out on +X with the Sun behind the camera.
+fn render(renderer: &mut Renderer) -> Frame {
+    let observer = Pose::new(Vector3::new(4.2e7, 0.0, 0.0), UnitQuaternion::identity());
+    let k = renderer.camera().intrinsics();
+    let moon = observer.position + k.ray(&Vector2::new(k.cx + 6.0, k.cy + 4.0)) * 3.844e8;
+    let scene = Scene {
+        earth_attitude: UnitQuaternion::identity(),
+        sun: observer.position - Vector3::z() * AU,
+        moon,
+        moon_attitude: UnitQuaternion::identity(),
+        targets: Vec::new(),
+    };
+    renderer.submit(&scene, &observer).unwrap().wait().unwrap()
+}
+
+/// The example camera cut down to 160 x 120 pixels.
+fn small_camera(psf: Psf) -> CameraConfig {
+    CameraConfig {
+        width: 160,
+        height: 120,
+        psf,
+        ..CameraConfig::example()
+    }
+}
+
+#[test]
+fn the_renderer_blurs_finer_pixels_from_camera_c1() {
+    // At C1, rendering finer and blurring keep the Moon's centroid, which a projection
+    // off by any part of a rendered pixel would move, and its total light. Rendering
+    // finer lowers the total by 1.8e-4 at 2 x 2 and 2.3e-4 at 3 x 3: a full Moon dims
+    // toward its edge along a curve, so a pixel's center reads slightly more than the
+    // pixel's average. At C0 the PSF is ignored, and finer rendering refused.
+    let _devices = common::device_lock();
+    for adapter in &ADAPTERS {
+        let mut reference =
+            renderer(adapter, small_camera(Psf::None), CameraTier::C1, Some(1)).unwrap();
+        let reference = render(&mut reference).radiance.unwrap();
+        for (psf, supersampling) in [
+            (Psf::None, Some(2)),
+            (sharp(), Some(1)),
+            (sharp(), None),
+            (sharp(), Some(3)),
+        ] {
+            let mut r = renderer(
+                adapter,
+                small_camera(psf.clone()),
+                CameraTier::C1,
+                supersampling,
+            )
+            .unwrap();
+            let name = format!("{}, {psf:?}, {supersampling:?}", r.adapter_info().name);
+            let image = render(&mut r).radiance.unwrap();
+            assert_eq!(image.len(), 160 * 120);
+            for k in 0..3 {
+                let moments = |image: &[[f32; 4]]| {
+                    let (mut total, mut u, mut v) = (0.0, 0.0, 0.0);
+                    for (i, p) in image.iter().enumerate() {
+                        let l = p[k] as f64;
+                        total += l;
+                        u += l * (i % 160) as f64;
+                        v += l * (i / 160) as f64;
+                    }
+                    (total, u / total, v / total)
+                };
+                let ((ta, ua, va), (tb, ub, vb)) = (moments(&reference), moments(&image));
+                assert!(
+                    (tb / ta - 1.0).abs() < 3e-4,
+                    "{name}: channel {k}: {tb} vs {ta}"
+                );
+                assert!(
+                    (ub - ua).abs() < 0.01 && (vb - va).abs() < 0.01,
+                    "{name}: channel {k}: centroid ({ub}, {vb}) vs ({ua}, {va})"
+                );
+            }
+        }
+
+        let c0 = |psf| {
+            let mut r = renderer(adapter, small_camera(psf), CameraTier::C0, None).unwrap();
+            render(&mut r).dn
+        };
+        assert_eq!(c0(Psf::None), c0(sharp()));
+        assert!(renderer(adapter, small_camera(sharp()), CameraTier::C0, Some(2)).is_err());
+        assert!(renderer(adapter, small_camera(sharp()), CameraTier::C1, Some(5)).is_err());
+    }
+}

crates/resolved/tests/gpu/sensor.rs

@@ -16,7 +16,7 @@ use resolved::{CameraConfig, Pose, Scene, Spectrum};
 /// 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, 1);
+    let targets = FrameTargets::new(ctx, config.width, config.height, 3, 1, 1, 0);
     ctx.queue.write_texture(
         wgpu::TexelCopyTextureInfo {
             texture: &targets.radiance,
@@ -159,7 +159,7 @@ fn sunlit_earth_through_the_whole_chain() {
     );
     let sensor = SensorUniforms::new(&config, &radiometry);
     for ctx in common::contexts() {
-        let targets = FrameTargets::new(ctx, config.width, config.height, 3, 1);
+        let targets = FrameTargets::new(ctx, config.width, config.height, 3, 1, 1, 0);
         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).encode(ctx, &mut encoder, &targets, &sensor, None);