From 49e0500007369b0834c93daa4c2c3ea3f4ad8375 Mon Sep 17 00:00:00 2001 From: Connor Johnstone Date: Thu, 15 Oct 2026 19:28:30 -0400 Subject: [PATCH] Lens distortion --- crates/resolved-cli/src/main.rs | 2 +- crates/resolved-cli/src/scene_file.rs | 24 +- crates/resolved/src/camera.rs | 237 +++++++++++++- crates/resolved/src/fidelity.rs | 5 +- crates/resolved/src/gpu/blur.rs | 85 +++-- crates/resolved/src/gpu/mod.rs | 1 + crates/resolved/src/gpu/profile.rs | 4 +- crates/resolved/src/gpu/sensor.rs | 19 +- crates/resolved/src/gpu/targets.rs | 100 ++++-- crates/resolved/src/gpu/warp.rs | 124 +++++++ crates/resolved/src/lib.rs | 2 +- crates/resolved/src/render.rs | 84 +++-- crates/resolved/src/shaders/blur.wgsl | 64 +++- crates/resolved/src/shaders/distortion.wgsl | 44 +++ crates/resolved/src/shaders/sensor.wgsl | 15 +- crates/resolved/tests/gpu/bodies.rs | 9 +- crates/resolved/tests/gpu/distortion.rs | 340 ++++++++++++++++++++ crates/resolved/tests/gpu/main.rs | 1 + crates/resolved/tests/gpu/mesh.rs | 9 +- crates/resolved/tests/gpu/pixel_area.rs | 2 +- crates/resolved/tests/gpu/psf.rs | 11 +- crates/resolved/tests/gpu/sensor.rs | 10 +- docs/optics.md | 80 ++++- readme.md | 7 +- 24 files changed, 1146 insertions(+), 133 deletions(-) create mode 100644 crates/resolved/src/gpu/warp.rs create mode 100644 crates/resolved/src/shaders/distortion.wgsl create mode 100644 crates/resolved/tests/gpu/distortion.rs diff --git a/crates/resolved-cli/src/main.rs b/crates/resolved-cli/src/main.rs index e2a96a1..a147f13 100644 --- a/crates/resolved-cli/src/main.rs +++ b/crates/resolved-cli/src/main.rs @@ -57,7 +57,7 @@ struct RenderArgs { #[arg(long, default_value = "L0")] lighting: LightingTier, /// Camera fidelity tier: C0, or C1 (adds pixels that collect light over their area, - /// for graded edges and continuous thin parts, and lens blur). + /// for graded edges and continuous thin parts, lens distortion 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 diff --git a/crates/resolved-cli/src/scene_file.rs b/crates/resolved-cli/src/scene_file.rs index 270d42f..cdf2ca5 100644 --- a/crates/resolved-cli/src/scene_file.rs +++ b/crates/resolved-cli/src/scene_file.rs @@ -10,8 +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, Psf, Renderer, Scene, Spectrum, - Target, + BrownConrady, CameraConfig, LoadOptions, Material, Mesh, ModelId, Pose, Psf, Renderer, Scene, + Spectrum, Target, }; use serde::Deserialize; @@ -64,6 +64,8 @@ pub struct CameraFile { /// 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, + /// Lens distortion from camera C1, in OpenCV's order: [k1, k2, p1, p2, k3]. + pub distortion: Option<[f64; 5]>, } /// One sigma, or one per channel. @@ -157,6 +159,9 @@ impl SceneFile { if let Some(t) = c.transmission { camera.transmission = Spectrum::constant(t); } + if let Some([k1, k2, p1, p2, k3]) = c.distortion { + camera.distortion = Some(BrownConrady { k1, k2, k3, p1, p2 }); + } camera.psf = match &c.psf_sigma { None => camera.psf, Some(Sigma::All(s)) if *s == 0.0 => Psf::None, @@ -410,6 +415,21 @@ mod tests { toml::from_str(&format!("{MINIMAL}\n[camera]\n{line}\n")).unwrap(); assert_eq!(file.camera().psf, psf, "{line}"); } + + let file: SceneFile = toml::from_str(&format!( + "{MINIMAL}\n[camera]\ndistortion = [-0.1, 0.02, 1e-3, -2e-3, 0.005]\n" + )) + .unwrap(); + assert_eq!( + file.camera().distortion, + Some(BrownConrady { + k1: -0.1, + k2: 0.02, + k3: 0.005, + p1: 1e-3, + p2: -2e-3, + }) + ); } #[test] diff --git a/crates/resolved/src/camera.rs b/crates/resolved/src/camera.rs index ab4995f..1888190 100644 --- a/crates/resolved/src/camera.rs +++ b/crates/resolved/src/camera.rs @@ -5,7 +5,7 @@ //! matching the camera frame in [`crate::scene`]: +X toward increasing `u`, +Y toward //! increasing `v`, +Z along the boresight. -use nalgebra::{Matrix4, UnitQuaternion, Vector2, Vector3}; +use nalgebra::{Matrix2, Matrix4, UnitQuaternion, Vector2, Vector3}; use crate::scene::Pose; use crate::spectrum::Spectrum; @@ -107,6 +107,76 @@ impl Psf { } } +/// Brown-Conrady lens distortion, in OpenCV's convention: radial terms `k1`, `k2`, `k3` +/// and tangential terms `p1`, `p2`, acting on normalized coordinates `x = X / Z`, +/// `y = Y / Z` in the camera frame, with `r^2 = x^2 + y^2`: +/// +/// ```text +/// x_d = x (1 + k1 r^2 + k2 r^4 + k3 r^6) + 2 p1 x y + p2 (r^2 + 2 x^2) +/// y_d = y (1 + k1 r^2 + k2 r^4 + k3 r^6) + p1 (r^2 + 2 y^2) + 2 p2 x y +/// ``` +/// +/// The intrinsics then map `(x_d, y_d)` to pixels. Coefficients from an OpenCV +/// calibration, whose `distCoeffs` are ordered `(k1, k2, p1, p2, k3)`, apply directly. +#[derive(Clone, Copy, Debug, Default, PartialEq)] +pub struct BrownConrady { + /// Radial, r^2. Positive moves points outward (pincushion), negative inward (barrel). + pub k1: f64, + /// Radial, r^4. + pub k2: f64, + /// Radial, r^6. + pub k3: f64, + /// Tangential, from a lens element decentered along y. + pub p1: f64, + /// Tangential, from a lens element decentered along x. + pub p2: f64, +} + +impl BrownConrady { + /// Where normalized point `p` lands after distortion. + pub fn distort(&self, p: &Vector2) -> Vector2 { + let (x, y) = (p.x, p.y); + let r2 = x * x + y * y; + let radial = 1.0 + r2 * (self.k1 + r2 * (self.k2 + r2 * self.k3)); + Vector2::new( + x * radial + 2.0 * self.p1 * x * y + self.p2 * (r2 + 2.0 * x * x), + y * radial + self.p1 * (r2 + 2.0 * y * y) + 2.0 * self.p2 * x * y, + ) + } + + /// The derivative of [`BrownConrady::distort`] at `p`: rows x_d and y_d, columns x + /// and y. Its determinant is how much the distortion stretches area there. + pub fn jacobian(&self, p: &Vector2) -> Matrix2 { + let (x, y) = (p.x, p.y); + let r2 = x * x + y * y; + let radial = 1.0 + r2 * (self.k1 + r2 * (self.k2 + r2 * self.k3)); + // d(radial)/d(r^2). + let q = self.k1 + r2 * (2.0 * self.k2 + 3.0 * r2 * self.k3); + let cross = 2.0 * x * y * q + 2.0 * self.p1 * x + 2.0 * self.p2 * y; + Matrix2::new( + radial + 2.0 * x * x * q + 2.0 * self.p1 * y + 6.0 * self.p2 * x, + cross, + cross, + radial + 2.0 * y * y * q + 6.0 * self.p1 * y + 2.0 * self.p2 * x, + ) + } + + /// The normalized point that [`BrownConrady::distort`] moves to `d`, by Newton's + /// method from `d`, or `None` if it does not converge to 1e-14 (1e-10 px for the + /// example camera's focal length). + pub fn undistort(&self, d: &Vector2) -> Option> { + let mut p = *d; + for _ in 0..50 { + let step = self.jacobian(&p).try_inverse()? * (self.distort(&p) - d); + p -= step; + if step.norm() < 1e-14 { + return Some(p); + } + } + None + } +} + /// A pinhole camera with radiometric optics and a detector. #[derive(Clone, Debug, PartialEq)] pub struct CameraConfig { @@ -128,6 +198,9 @@ pub struct CameraConfig { pub transmission: Spectrum, /// How the optics blur the image. pub psf: Psf, + /// How the lens bends straight lines, from camera tier C1; `None` for a perfect + /// pinhole. + pub distortion: Option, /// The detector's color channels. pub channels: Vec, /// Electrons to DN. @@ -152,6 +225,7 @@ impl CameraConfig { /// - 20,000 electron full well, 12-bit ADC, 5 electrons per DN, 64 DN bias. /// - 0.1 ms exposure: a white Lambertian surface facing the Sun fills about half the /// well (48-57 % across the channels), so Earth and the target are both unsaturated. + /// - No lens distortion. /// - Mounted with no offset or rotation on the observer. pub fn example() -> Self { let channel = |name: &str, center: f64| Channel { @@ -170,6 +244,7 @@ impl CameraConfig { psf: Psf::Gaussian { sigma: vec![0.5; 3], }, + distortion: None, channels: vec![ channel("r", 620e-9), channel("g", 540e-9), @@ -215,6 +290,15 @@ impl CameraConfig { return fail("PSF sigma must be positive and at most 15 pixels"); } } + if let Some(distortion) = &self.distortion { + let b = distortion; + if ![b.k1, b.k2, b.k3, b.p1, b.p2].iter().all(|x| x.is_finite()) { + return fail("distortion coefficients must be finite"); + } + if !self.intrinsics().distortion_invertible() { + return fail("the lens distortion folds the image over itself"); + } + } 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"); @@ -225,7 +309,7 @@ impl CameraConfig { Ok(()) } - /// Pinhole projection parameters in pixels. + /// Projection parameters in pixels, with the lens distortion. pub fn intrinsics(&self) -> Intrinsics { let f = self.focal_length / self.pixel_pitch; let c = self @@ -238,11 +322,13 @@ impl CameraConfig { cy: c.y, width: self.width, height: self.height, + distortion: self.distortion, } } } -/// Pinhole projection in pixels: `u = fx x / z + cx`, `v = fy y / z + cy`. +/// Projection in pixels: normalized coordinates `x / z`, `y / z`, distorted by +/// `distortion` if any, then `u = fx x_d + cx`, `v = fy y_d + cy`. #[derive(Clone, Copy, Debug, PartialEq)] pub struct Intrinsics { /// Focal length in pixels along `u`. @@ -257,22 +343,71 @@ pub struct Intrinsics { pub width: u32, /// Image height in pixels. pub height: u32, + /// The lens distortion; `None` for a perfect pinhole. + pub distortion: Option, } impl Intrinsics { /// Where a camera-frame point lands in the image, or `None` if it is not in front of /// the camera. The result may lie outside the image. pub fn project(&self, p: &Vector3) -> Option> { - (p.z > 0.0) - .then(|| Vector2::new(self.fx * p.x / p.z + self.cx, self.fy * p.y / p.z + self.cy)) + (p.z > 0.0).then(|| { + let mut n = Vector2::new(p.x / p.z, p.y / p.z); + if let Some(d) = &self.distortion { + n = d.distort(&n); + } + self.from_normalized(&n) + }) } /// Unit direction in the camera frame of the ray through pixel coordinates `uv`. + /// + /// With distortion, this inverts it by Newton's method, to 1e-10 px for the example + /// camera. Where the distortion cannot be inverted, which + /// [`CameraConfig::validate`] rules out within the image, the result is the + /// method's last estimate. pub fn ray(&self, uv: &Vector2) -> Vector3 { - Vector3::new((uv.x - self.cx) / self.fx, (uv.y - self.cy) / self.fy, 1.0).normalize() + let n = self.to_normalized(uv); + let n = match &self.distortion { + Some(d) => d.undistort(&n).unwrap_or(n), + None => n, + }; + Vector3::new(n.x, n.y, 1.0).normalize() } - /// Camera frame to wgpu clip space, with reversed-Z depth and no far plane. + /// Pixel coordinates of normalized (distorted) coordinates `n`. + pub fn from_normalized(&self, n: &Vector2) -> Vector2 { + Vector2::new(self.fx * n.x + self.cx, self.fy * n.y + self.cy) + } + + /// Normalized (distorted) coordinates of pixel coordinates `uv`. + pub fn to_normalized(&self, uv: &Vector2) -> Vector2 { + Vector2::new((uv.x - self.cx) / self.fx, (uv.y - self.cy) / self.fy) + } + + /// Whether the distortion can be inverted across the image: at every pixel along + /// its border and on a grid 16 pixels apart inside, it converges and leaves area + /// positive (it does not fold the image over). + pub fn distortion_invertible(&self) -> bool { + let Some(d) = &self.distortion else { + return true; + }; + let (w, h) = (self.width as f64, self.height as f64); + let along = |n: f64, step: f64| { + (0..=(n / step).ceil() as u32).map(move |i| (i as f64 * step).min(n)) + }; + let border = along(w, 1.0) + .flat_map(|u| [Vector2::new(u, 0.0), Vector2::new(u, h)]) + .chain(along(h, 1.0).flat_map(|v| [Vector2::new(0.0, v), Vector2::new(w, v)])); + let grid = along(w, 16.0).flat_map(|u| along(h, 16.0).map(move |v| Vector2::new(u, v))); + border.chain(grid).all(|uv| { + d.undistort(&self.to_normalized(&uv)) + .is_some_and(|p| d.jacobian(&p).determinant() > 0.0) + }) + } + + /// Camera frame to wgpu clip space, with reversed-Z depth and no far plane: the + /// pinhole part of the projection, without distortion, which a rasterizer can apply. /// /// Depth is `near / z`: 1 at `near`, falling toward 0 at infinity, which spends f32 /// precision evenly in log distance from a 10 m target out to the Sun. Clip `x` and @@ -397,6 +532,94 @@ mod tests { } } + /// Typical coefficients for a small lens, with all five terms. + fn distorted() -> CameraConfig { + CameraConfig { + principal_point: Some(Vector2::new(955.25, 543.5)), + distortion: Some(BrownConrady { + k1: -0.12, + k2: 0.05, + k3: 0.02, + p1: 8e-4, + p2: -5e-4, + }), + ..CameraConfig::example() + } + } + + #[test] + fn distortion_matches_opencv() { + // cv2.projectPoints (OpenCV 5.0.0) with distCoeffs (k1, k2, p1, p2, k3). + let k = distorted().intrinsics(); + for (p, uv) in [ + ([0.3, -0.2, 10.0], [1227.9119321483472, 361.7308937192836]), + ( + [-1.0e5, 4.0e5, 7.0e6], + [825.4043611695719, 1062.8447074552932], + ), + ([-0.9, -0.5, 8.0], [-65.54679093052266, -23.44704925433075]), + ([2.0, 1.0, 15.0], [2164.0357505092074, 1148.104996466725]), + ] { + let got = k.project(&Vector3::from(p)).unwrap(); + assert_relative_eq!(got, Vector2::from(uv), epsilon = 1e-9); + } + } + + #[test] + fn distorted_ray_inverts_project() { + // Across the image and a little past it, to 1e-9 px. + let k = distorted().intrinsics(); + for i in 0..=20 { + for j in 0..=20 { + let uv = Vector2::new(-40.0 + 100.0 * i as f64, -30.0 + 57.0 * j as f64); + let back = k.project(&(k.ray(&uv) * 12.0)).unwrap(); + assert!((back - uv).norm() < 1e-9, "{uv}: {back}"); + } + } + } + + #[test] + fn distortion_jacobian_matches_differences() { + let d = distorted().distortion.unwrap(); + let p = Vector2::new(0.08, -0.05); + let h = 1e-7; + let j = d.jacobian(&p); + for (col, e) in [Vector2::x(), Vector2::y()].into_iter().enumerate() { + let diff = (d.distort(&(p + e * h)) - d.distort(&(p - e * h))) / (2.0 * h); + assert_relative_eq!(j.column(col).into_owned(), diff, epsilon = 1e-8); + } + } + + #[test] + fn radial_distortion_moves_corners_outward() { + // k1 = 0.1 on the example camera: the corner, 0.0146 in r^2, moves out by + // k1 r^2 of its 1101 px from the center, 1.6 px. + let mut camera = CameraConfig::example(); + camera.distortion = Some(BrownConrady { + k1: 0.1, + ..Default::default() + }); + let k = camera.intrinsics(); + let pinhole = CameraConfig::example().intrinsics(); + let corner = Vector3::new(960.0 / k.fx, 540.0 / k.fy, 1.0); + let moved = k.project(&corner).unwrap() - pinhole.project(&corner).unwrap(); + let r2 = corner.x * corner.x + corner.y * corner.y; + assert_relative_eq!(moved.norm(), 0.1 * r2 * 1101.4, max_relative = 1e-3); + assert_relative_eq!(moved.norm(), 1.61, epsilon = 0.01); + } + + #[test] + fn validate_rejects_folding_distortion() { + distorted().validate().unwrap(); + let mut camera = CameraConfig::example(); + // 1 + 3 k1 r^2 turns negative before the corners. + camera.distortion = Some(BrownConrady { + k1: -30.0, + ..Default::default() + }); + assert!(camera.validate().is_err()); + } + #[test] fn clip_matrix_agrees_with_projection() { let mut camera = CameraConfig::example(); diff --git a/crates/resolved/src/fidelity.rs b/crates/resolved/src/fidelity.rs index 396ecad..b89c3b5 100644 --- a/crates/resolved/src/fidelity.rs +++ b/crates/resolved/src/fidelity.rs @@ -39,9 +39,10 @@ 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, and the lens's blur (the camera's + /// features thinner than a pixel stay continuous, the lens's distortion (the camera's + /// [`BrownConrady`](crate::BrownConrady)) and its 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 + /// default, distorted, 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. diff --git a/crates/resolved/src/gpu/blur.rs b/crates/resolved/src/gpu/blur.rs index e72b418..85dbd0a 100644 --- a/crates/resolved/src/gpu/blur.rs +++ b/crates/resolved/src/gpu/blur.rs @@ -1,11 +1,17 @@ -//! The pass that blurs the rendered radiance by the camera's point spread function and -//! averages it into the image's pixels. +//! The pass that makes the focal plane's image: it applies lens distortion, 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. //! +//! With lens distortion the render is undistorted, and the rows pass reads it through +//! the distortion as it loads its tile (see [`super::warp`]), rather than a separate +//! pass writing the distorted image out and this one reading it back: at 2 x 2 +//! rendering that measured 0.7 to 0.8 ms per 1920 x 1080 frame on an RTX 3060. +//! //! 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 @@ -21,6 +27,7 @@ use super::profile::PassTimestamps; use super::targets::FrameTargets; +use super::warp::WarpUniforms; use super::GpuContext; use crate::camera::Psf; @@ -41,27 +48,40 @@ struct Params { scale: u32, radius: u32, margin: u32, - unused: u32, + warp: u32, + focal: [u32; 2], + unused: [u32; 2], } -/// The compiled rows and columns pipelines and the weights. +/// The compiled rows and columns pipelines, the weights, and each pass's uniforms. #[derive(Debug)] pub struct BlurPass { line: u32, rows: wgpu::ComputePipeline, columns: wgpu::ComputePipeline, weights: wgpu::Buffer, - params: wgpu::Buffer, + /// [`Params`] for the rows pass, then the columns pass. + params: [wgpu::Buffer; 2], + /// [`WarpUniforms`], zero without distortion. + warp: 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 { + /// over radiance rendered `scale` times finer than the image of `size` pixels, with + /// a margin of `psf.radius(scale)`. With `warp`, the render is undistorted and the + /// rows pass reads the focal plane through the distortion. `None` with no PSF, a + /// `scale` of 1 and no distortion, where there is nothing to do. + pub fn new( + ctx: &GpuContext, + psf: &Psf, + channels: usize, + scale: u32, + size: [u32; 2], + warp: Option<&WarpUniforms>, + ) -> Option { use wgpu::util::DeviceExt; - if *psf == Psf::None && scale == 1 { + if *psf == Psf::None && scale == 1 && warp.is_none() { return None; } let per_channel: Vec> = (0..channels) @@ -71,17 +91,24 @@ impl BlurPass { .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 { + let params = |warp: bool| Params { scale, radius, margin: radius, - unused: 0, + warp: warp.into(), + focal: size.map(|n| n * scale + 2 * radius), + unused: [0; 2], }; + let source = [ + include_str!("../shaders/distortion.wgsl"), + include_str!("../shaders/blur.wgsl"), + ] + .concat(); let module = ctx .device .create_shader_module(wgpu::ShaderModuleDescriptor { label: Some("blur"), - source: wgpu::ShaderSource::Wgsl(include_str!("../shaders/blur.wgsl").into()), + source: wgpu::ShaderSource::Wgsl(source.into()), }); let pipeline = |entry| { ctx.device @@ -114,9 +141,16 @@ impl BlurPass { bytemuck::cast_slice(&weights), wgpu::BufferUsages::STORAGE, ), - params: buffer( - "blur params", - bytemuck::bytes_of(¶ms), + params: [params(warp.is_some()), params(false)].map(|p| { + buffer( + "blur params", + bytemuck::bytes_of(&p), + wgpu::BufferUsages::UNIFORM, + ) + }), + warp: buffer( + "blur warp", + bytemuck::bytes_of(&warp.copied().unwrap_or_default()), wgpu::BufferUsages::UNIFORM, ), }) @@ -134,19 +168,26 @@ impl BlurPass { let [rows, image] = targets .blurred .as_ref() - .expect("targets allocated with a margin or a scale"); + .expect("targets allocated with a margin, a scale or distortion"); 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)), + for (pipeline, source, destination, groups, params) in [ + ( + &self.rows, + &targets.radiance, + rows, + (w.div_ceil(line), h), + &self.params[0], + ), ( &self.columns, rows, image, (w, image.height().div_ceil(line)), + &self.params[1], ), ] { let bind_group = ctx.device.create_bind_group(&wgpu::BindGroupDescriptor { @@ -171,7 +212,11 @@ impl BlurPass { }, wgpu::BindGroupEntry { binding: 3, - resource: self.params.as_entire_binding(), + resource: params.as_entire_binding(), + }, + wgpu::BindGroupEntry { + binding: 4, + resource: self.warp.as_entire_binding(), }, ], }); diff --git a/crates/resolved/src/gpu/mod.rs b/crates/resolved/src/gpu/mod.rs index d22a320..1632575 100644 --- a/crates/resolved/src/gpu/mod.rs +++ b/crates/resolved/src/gpu/mod.rs @@ -9,5 +9,6 @@ pub mod profile; pub mod sensor; pub mod shadow; pub mod targets; +pub mod warp; pub use context::{adapters, AdapterChoice, GpuContext}; diff --git a/crates/resolved/src/gpu/profile.rs b/crates/resolved/src/gpu/profile.rs index 2ceac64..f1c3f8e 100644 --- a/crates/resolved/src/gpu/profile.rs +++ b/crates/resolved/src/gpu/profile.rs @@ -71,8 +71,8 @@ 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). + /// Lens distortion, the point spread function and averaging finer rendering into + /// pixels (zero work below camera C1). pub blur: f64, /// Radiance to DN. pub sensor: f64, diff --git a/crates/resolved/src/gpu/sensor.rs b/crates/resolved/src/gpu/sensor.rs index b28f0d4..b035367 100644 --- a/crates/resolved/src/gpu/sensor.rs +++ b/crates/resolved/src/gpu/sensor.rs @@ -6,6 +6,7 @@ use super::profile::PassTimestamps; use super::targets::FrameTargets; +use super::warp::DistortionUniform; use super::GpuContext; use crate::camera::CameraConfig; use crate::radiometry::Radiometry; @@ -17,12 +18,14 @@ pub struct SensorUniforms { intrinsics: [f32; 4], detector: [f32; 4], dims: [u32; 4], + distortion: DistortionUniform, } 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 { + /// The constants for `camera`, whose channels `radiometry` was built from, with its + /// lens distortion if `distorted` (camera C1 and up). They do not change from frame + /// to frame. + pub fn new(camera: &CameraConfig, radiometry: &Radiometry, distorted: bool) -> Self { let k = camera.intrinsics(); let d = &camera.detector; Self { @@ -40,6 +43,7 @@ impl SensorUniforms { camera.channels.len() as u32, d.dn_max(), ], + distortion: DistortionUniform::new(k.distortion.as_ref().filter(|_| distorted)), } } } @@ -60,7 +64,14 @@ impl SensorPass { .device .create_shader_module(wgpu::ShaderModuleDescriptor { label: Some("sensor"), - source: wgpu::ShaderSource::Wgsl(include_str!("../shaders/sensor.wgsl").into()), + source: wgpu::ShaderSource::Wgsl( + [ + include_str!("../shaders/distortion.wgsl"), + include_str!("../shaders/sensor.wgsl"), + ] + .concat() + .into(), + ), }); let pipeline = ctx .device diff --git a/crates/resolved/src/gpu/targets.rs b/crates/resolved/src/gpu/targets.rs index 39d2eaf..8f8dfc5 100644 --- a/crates/resolved/src/gpu/targets.rs +++ b/crates/resolved/src/gpu/targets.rs @@ -12,6 +12,52 @@ pub const DEPTH_FORMAT: wgpu::TextureFormat = wgpu::TextureFormat::Depth32Float; /// Depth test for reversed Z. pub const DEPTH_COMPARE: wgpu::CompareFunction = wgpu::CompareFunction::Greater; +/// The sizes of one frame's images, for [`FrameTargets::new`]. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub struct TargetLayout { + /// Image width in pixels. + pub width: u32, + /// Image height in pixels. + pub height: u32, + /// Number of camera channels. + pub channels: u32, + /// Mesh samples per rendered pixel: 1, 2 or 4. + pub samples: u32, + /// Rendered pixels along each axis per image pixel. + pub scale: u32, + /// Rendered pixels past the image on each side, for the blur. + pub margin: u32, + /// With lens distortion, the size of the undistorted render that the focal plane is + /// read from; `None` without, where the render is the focal plane. + pub undistorted: Option<[u32; 2]>, +} + +impl TargetLayout { + /// A `width` x `height` image of `channels` channels, rendered one pixel per pixel + /// with one mesh sample, no margin and no distortion. + pub fn new(width: u32, height: u32, channels: u32) -> Self { + Self { + width, + height, + channels, + samples: 1, + scale: 1, + margin: 0, + undistorted: None, + } + } + + /// The focal plane's size in rendered pixels: the image at `scale`, with the margin. + pub fn focal_plane(&self) -> [u32; 2] { + [self.width, self.height].map(|n| n * self.scale + 2 * self.margin) + } + + /// The size the meshes and bodies render at. + pub fn rendered(&self) -> [u32; 2] { + self.undistorted.unwrap_or_else(|| self.focal_plane()) + } +} + /// The images one frame renders into, sized for one camera. /// /// Target meshes are drawn first, into their own radiance and depth with @@ -20,8 +66,10 @@ pub const DEPTH_COMPARE: wgpu::CompareFunction = wgpu::CompareFunction::Greater; /// 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 +/// light from just outside the frame can blur into it. With lens distortion they are +/// rendered undistorted, on the area [`TargetLayout::undistorted`] gives. The blur +/// passes then read the focal plane from them, through the distortion if any, blur it +/// and average it into the image's pixels, in `blurred`. The sensor pass reads /// [`FrameTargets::sensed`]. #[derive(Debug)] pub struct FrameTargets { @@ -37,7 +85,7 @@ pub struct FrameTargets { /// several points. pub samples: u32, /// The meshes' channel radiance, [`RADIANCE_FORMAT`], `samples` per pixel. This and - /// the next three are rendered at `scale` and include the margin. + /// the next three are the size [`TargetLayout::rendered`] gives. pub mesh_radiance: wgpu::Texture, /// The meshes' depth, [`DEPTH_FORMAT`], `samples` per pixel. pub mesh_depth: wgpu::Texture, @@ -46,9 +94,9 @@ 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`]. + /// With a blur, a `scale` above 1 or distortion: the focal plane 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, @@ -58,20 +106,20 @@ pub struct FrameTargets { } impl FrameTargets { - /// 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); + /// Allocates the images `layout` describes. The blur's images are allocated with a + /// margin, a scale above 1 or distortion. + pub fn new(ctx: &GpuContext, layout: &TargetLayout) -> Self { + let TargetLayout { + width, + height, + channels, + samples, + scale, + margin, + .. + } = *layout; + let [outer_width, outer_height] = layout.rendered(); + let focal = layout.focal_plane(); let sized = |label, format, usage, sample_count, width, height| { ctx.device.create_texture(&wgpu::TextureDescriptor { label: Some(label), @@ -103,7 +151,7 @@ impl FrameTargets { } else { wgpu::TextureUsages::empty() }; - let blur = |label, height| { + let storage = |label, width, height| { sized( label, RADIANCE_FORMAT, @@ -131,8 +179,12 @@ impl FrameTargets { 1, ), depth: texture("depth", DEPTH_FORMAT, usage, 1), - blurred: (margin > 0 || scale > 1) - .then(|| [blur("blurred rows", outer_height), blur("blurred", height)]), + blurred: (margin > 0 || scale > 1 || layout.undistorted.is_some()).then(|| { + [ + storage("blurred rows", width, focal[1]), + storage("blurred", width, height), + ] + }), channels, dn: ctx.device.create_buffer(&wgpu::BufferDescriptor { label: Some("dn"), @@ -165,7 +217,7 @@ impl FrameTargets { } /// The radiance image the sensor pass reads, `width` x `height`: the blurred image, - /// or without a blur, `radiance`. + /// or without one, `radiance`. pub fn sensed(&self) -> &wgpu::Texture { match &self.blurred { Some([_, image]) => image, diff --git a/crates/resolved/src/gpu/warp.rs b/crates/resolved/src/gpu/warp.rs new file mode 100644 index 0000000..7499507 --- /dev/null +++ b/crates/resolved/src/gpu/warp.rs @@ -0,0 +1,124 @@ +//! Lens distortion: where the undistorted render lies, and what the focal-plane pass +//! needs to read it through the distortion. +//! +//! Meshes are rasterized with a linear projection, which keeps straight edges straight, +//! so the frame is drawn without distortion and then warped once, by the blur's rows +//! pass as it reads the render ([`super::blur`]). For each focal-plane pixel the shader +//! inverts the distortion by Newton's method, which costs a few multiplications and +//! saves storing a map of 8 bytes per rendered pixel, and reads the render there +//! bilinearly. The render is sized and placed by [`WarpLayout::new`] to cover every +//! point the focal plane reads. + +use nalgebra::Vector2; + +use crate::camera::{BrownConrady, Intrinsics}; +use crate::Error; + +/// `Distortion` in `distortion.wgsl`. +#[repr(C)] +#[derive(Clone, Copy, Debug, Default, PartialEq, bytemuck::Pod, bytemuck::Zeroable)] +pub struct DistortionUniform { + k: [f32; 4], + p: [f32; 4], +} + +impl DistortionUniform { + /// `distortion`, or a pinhole for `None`. + pub fn new(distortion: Option<&BrownConrady>) -> Self { + match distortion { + None => Self::default(), + Some(d) => Self { + k: [d.k1 as f32, d.k2 as f32, d.k3 as f32, 1.0], + p: [d.p1 as f32, d.p2 as f32, 0.0, 0.0], + }, + } + } +} + +/// What the focal-plane pass needs to read the render through the distortion; the +/// layout matches `Warp` in `blur.wgsl`. +#[repr(C)] +#[derive(Clone, Copy, Debug, Default, PartialEq, bytemuck::Pod, bytemuck::Zeroable)] +pub struct WarpUniforms { + distortion: DistortionUniform, + image: [f32; 4], + offset: [f32; 4], +} + +impl WarpUniforms { + /// For an image with intrinsics `k`, whose distortion it applies, and an undistorted + /// render placed by `layout`. + pub fn new(k: &Intrinsics, layout: &WarpLayout) -> Self { + Self { + distortion: DistortionUniform::new(k.distortion.as_ref()), + image: [k.fx, k.fy, k.cx, k.cy].map(|x| x as f32), + offset: [layout.offset[0] as f32, layout.offset[1] as f32, 0.0, 0.0], + } + } +} + +/// Where the undistorted render lies relative to the focal plane. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub struct WarpLayout { + /// The render's size, in rendered pixels. + pub size: [u32; 2], + /// Whole rendered pixels from a focal-plane pixel to the render pixel that sees the + /// same point when there is no distortion. + pub offset: [u32; 2], +} + +impl WarpLayout { + /// The render for an image with intrinsics `k` (and its distortion), rendered + /// `scale` times finer with a focal-plane margin of `margin` rendered pixels: large + /// enough that every focal-plane pixel reads inside it, with a pixel to spare for + /// the bilinear read. Errors if the distortion cannot be inverted out to the margin. + pub fn new(k: &Intrinsics, scale: u32, margin: u32) -> Result { + let d = k.distortion.unwrap_or_default(); + let s = scale as f64; + // The focal plane's extent in image pixels. + let pad = margin as f64 / s; + let (u0, u1) = (-pad, k.width as f64 + pad); + let (v0, v1) = (-pad, k.height as f64 + pad); + let along = |a: f64, b: f64, step: f64| { + let n = ((b - a) / step).ceil() as u32; + (0..=n).map(move |i| (a + i as f64 * step).min(b)) + }; + let border = along(u0, u1, 1.0) + .flat_map(|u| [Vector2::new(u, v0), Vector2::new(u, v1)]) + .chain(along(v0, v1, 1.0).flat_map(|v| [Vector2::new(u0, v), Vector2::new(u1, v)])); + let grid = + along(u0, u1, 16.0).flat_map(|u| along(v0, v1, 16.0).map(move |v| Vector2::new(u, v))); + let (mut low, mut high) = (Vector2::zeros(), Vector2::::zeros()); + for uv in border.chain(grid) { + let distorted = k.to_normalized(&uv); + let n = d.undistort(&distorted).ok_or_else(|| { + Error::InvalidCamera(format!( + "the lens distortion cannot be inverted at pixel ({}, {})", + uv.x, uv.y + )) + })?; + // In rendered pixels. + let shift = (n - distorted).component_mul(&Vector2::new(k.fx, k.fy)) * s; + low = low.inf(&shift); + high = high.sup(&shift); + } + let focal = [k.width, k.height].map(|n| n * scale + 2 * margin); + let offset = [0, 1].map(|i| (-low[i]).ceil() as u32 + 1); + let size = [0, 1].map(|i| focal[i] + offset[i] + high[i].ceil() as u32 + 1); + Ok(Self { size, offset }) + } + + /// The pinhole intrinsics of the undistorted render, for image intrinsics `k`. + pub fn intrinsics(&self, k: &Intrinsics, scale: u32, margin: u32) -> Intrinsics { + let (s, m) = (scale as f64, margin as f64); + Intrinsics { + fx: k.fx * s, + fy: k.fy * s, + cx: k.cx * s + m + self.offset[0] as f64, + cy: k.cy * s + m + self.offset[1] as f64, + width: self.size[0], + height: self.size[1], + distortion: None, + } + } +} diff --git a/crates/resolved/src/lib.rs b/crates/resolved/src/lib.rs index 212f893..596141f 100644 --- a/crates/resolved/src/lib.rs +++ b/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, Psf}; +pub use camera::{BrownConrady, CameraConfig, Channel, Detector, Intrinsics, Psf}; pub use earth::EarthTexture; pub use error::Error; pub use fidelity::{CameraTier, Fidelity, LightingTier}; diff --git a/crates/resolved/src/render.rs b/crates/resolved/src/render.rs index c805afb..9f7d74d 100644 --- a/crates/resolved/src/render.rs +++ b/crates/resolved/src/render.rs @@ -25,7 +25,8 @@ use crate::gpu::mesh::{GpuMesh, MeshDraw, MeshPass, ShadowDraw}; use crate::gpu::profile::{PassTimings, Profiler}; use crate::gpu::sensor::{SensorPass, SensorUniforms}; use crate::gpu::shadow::{CasterDraw, ShadowPass}; -use crate::gpu::targets::{FrameTargets, DEPTH_FORMAT, RADIANCE_FORMAT}; +use crate::gpu::targets::{FrameTargets, TargetLayout, DEPTH_FORMAT, RADIANCE_FORMAT}; +use crate::gpu::warp::{WarpLayout, WarpUniforms}; use crate::gpu::{AdapterChoice, GpuContext}; use crate::model::{Mesh, ModelId}; use crate::moon::MoonTexture; @@ -114,7 +115,13 @@ pub struct Renderer { targets: FrameTargets, bodies: BodiesPass, meshes: MeshPass, + /// Whether the frame is distorted (camera C1 and a lens with distortion). + distorted: bool, blur: Option, + /// The pinhole projection the meshes and bodies render with: finer than the image + /// by the scale, with the margin, and with distortion, undistorted and placed for + /// the warp. + render_intrinsics: Intrinsics, shadows: ShadowPass, earthshine: EarthshinePass, sensor: SensorPass, @@ -139,7 +146,7 @@ impl Renderer { body_radiometry.bodies_eclipsed = options.fidelity.lighting >= LightingTier::L1; // 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. + // PSF reaches, then distorted, blurred and averaged into pixels. let area = options.fidelity.camera >= CameraTier::C1; let scale = match (options.fidelity.camera, options.supersampling) { (CameraTier::C0, None | Some(1)) => 1, @@ -169,9 +176,47 @@ impl Renderer { (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); + // From camera C1, lens distortion: render undistorted, and read the focal plane + // through the distortion. + let image = camera.intrinsics(); + let warp_layout = match &image.distortion { + Some(_) if area => Some(WarpLayout::new(&image, scale, margin)?), + _ => None, + }; + let warp = warp_layout.map(|l| WarpUniforms::new(&image, &l)); + let blur = BlurPass::new( + &ctx, + psf, + camera.channels.len(), + scale, + [camera.width, camera.height], + warp.as_ref(), + ); + let layout = TargetLayout { + samples, + scale, + margin, + undistorted: warp_layout.map(|l| l.size), + ..TargetLayout::new(camera.width, camera.height, camera.channels.len() as u32) + }; + let rendered = layout.rendered(); + // The pinhole projection the meshes and bodies render with. + let render_intrinsics = match &warp_layout { + Some(l) => l.intrinsics(&image, scale, margin), + None => { + let (f, m) = (scale as f64, margin as f64); + Intrinsics { + fx: image.fx * f, + fy: image.fy * f, + cx: image.cx * f + m, + cy: image.cy * f + m, + width: rendered[0], + height: rendered[1], + distortion: None, + } + } + }; let largest = ctx.device.limits().max_texture_dimension_2d; if rendered.iter().any(|&n| n > largest) { return Err(Error::Unsupported(format!( @@ -190,15 +235,7 @@ 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, - ); + let targets = FrameTargets::new(&ctx, &layout); 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}", @@ -260,7 +297,9 @@ impl Renderer { shadows: ShadowPass::new(&ctx, options.shadow_map_size.unwrap_or(4096)), earthshine: EarthshinePass::new(&ctx), sensor: SensorPass::new(&ctx), - sensor_uniforms: SensorUniforms::new(&camera, &radiometry), + sensor_uniforms: SensorUniforms::new(&camera, &radiometry, warp.is_some()), + distorted: warp.is_some(), + render_intrinsics, ctx, fidelity: options.fidelity, camera, @@ -347,7 +386,7 @@ impl Renderer { } self.radiometry.electrons_per_radiance *= seconds / self.camera.exposure; self.camera.exposure = seconds; - self.sensor_uniforms = SensorUniforms::new(&self.camera, &self.radiometry); + self.sensor_uniforms = SensorUniforms::new(&self.camera, &self.radiometry, self.distorted); Ok(()) } @@ -385,20 +424,7 @@ impl Renderer { ))); } let camera = observer.then(&self.camera.mount); - // 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 intrinsics = self.render_intrinsics; let near = self.camera.near_clip; let occulters = [ (Vector3::zeros(), crate::bodies::EARTH_EQUATORIAL_RADIUS), diff --git a/crates/resolved/src/shaders/blur.wgsl b/crates/resolved/src/shaders/blur.wgsl index fec0a17..20f9106 100644 --- a/crates/resolved/src/shaders/blur.wgsl +++ b/crates/resolved/src/shaders/blur.wgsl @@ -1,11 +1,19 @@ -// 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. +// The focal plane: lens distortion, the point spread function and the pixels' area. +// The rendered radiance, `scale` times finer than the image, is 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. +// (blur.rs). The focal plane extends `margin` rendered pixels past the image on both +// sides. The rows pass leaves the columns as rendered, margin included, for the +// columns pass. +// +// With lens distortion the render is undistorted, and the rows pass reads the focal +// plane from it through the distortion as it fills its tile: for each focal-plane +// pixel, the undistorted point the lens bends onto its center, read bilinearly. +// Radiance is unchanged along a ray, so moving it is all distortion does to an +// extended scene; the sensor accounts for how it stretches area. The code that +// includes this file first includes distortion.wgsl. struct Params { // Rendered pixels per image pixel along each axis. @@ -14,7 +22,20 @@ struct Params { radius: u32, // Rendered pixels past the image on each side, at least the radius. margin: u32, - unused: u32, + // 1 if this pass reads the focal plane through the distortion. + warp: u32, + // The focal plane's size in rendered pixels. + focal: vec2, + unused: vec2, +} + +struct Warp { + distortion: Distortion, + // fx, fy, cx, cy of the image, in image pixels. + image: vec4, + // xy: whole rendered pixels from a focal-plane pixel to the render pixel that sees + // the same point without distortion. + offset: vec4, } // Per tap, the weight of each channel. @@ -22,6 +43,30 @@ struct Params { @group(0) @binding(1) var source: texture_2d; @group(0) @binding(2) var destination: texture_storage_2d; @group(0) @binding(3) var p: Params; +@group(0) @binding(4) var w: Warp; + +fn texel(at: vec2) -> vec4 { + let size = vec2(textureDimensions(source)); + return textureLoad(source, clamp(at, vec2(0), size - 1), 0); +} + +// The focal plane's radiance at pixel `at`, read through the distortion from the +// undistorted render. +fn warped(at: vec2) -> vec4 { + let scale = f32(p.scale); + // This pixel's center in image pixels, then normalized and distorted. + let uv = (vec2(at) + 0.5 - f32(p.margin)) / scale; + let distorted = (uv - w.image.zw) / w.image.xy; + let shift = (undistort(w.distortion, distorted) - distorted) * w.image.xy * scale; + // In the render's texels, whose centers are at whole numbers. + let x = vec2(at) + w.offset.xy + shift; + let base = floor(x); + let f = x - base; + let i = vec2(base); + let top = mix(texel(i), texel(i + vec2(1, 0)), f.x); + let bottom = mix(texel(i + vec2(0, 1)), texel(i + vec2(1, 1)), f.x); + return mix(top, bottom, f.y); +} // 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 @@ -41,8 +86,11 @@ fn blur(start: vec2, i: u32, axis: vec2) { 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); + if p.warp == 1u { + tile[j] = warped(min(first + j * axis, p.focal - 1u)); + } else { + tile[j] = textureLoad(source, min(first + j * axis, source_size - 1u), 0); + } } workgroupBarrier(); let at = start + i * axis; diff --git a/crates/resolved/src/shaders/distortion.wgsl b/crates/resolved/src/shaders/distortion.wgsl new file mode 100644 index 0000000..1ff8d43 --- /dev/null +++ b/crates/resolved/src/shaders/distortion.wgsl @@ -0,0 +1,44 @@ +// Brown-Conrady lens distortion (BrownConrady in camera.rs), in f32, on normalized +// coordinates x = X / Z, y = Y / Z. + +struct Distortion { + // k1, k2, k3, and 1 if the lens distorts at all, 0 for a pinhole. + k: vec4, + // p1, p2. + p: vec4, +} + +fn distort(d: Distortion, n: vec2) -> vec2 { + let r2 = dot(n, n); + let radial = 1.0 + r2 * (d.k.x + r2 * (d.k.y + r2 * d.k.z)); + return vec2( + n.x * radial + 2.0 * d.p.x * n.x * n.y + d.p.y * (r2 + 2.0 * n.x * n.x), + n.y * radial + d.p.x * (r2 + 2.0 * n.y * n.y) + 2.0 * d.p.y * n.x * n.y, + ); +} + +// The derivative of distort at n; columns x and y. +fn distortion_jacobian(d: Distortion, n: vec2) -> mat2x2 { + let r2 = dot(n, n); + let radial = 1.0 + r2 * (d.k.x + r2 * (d.k.y + r2 * d.k.z)); + let q = d.k.x + r2 * (2.0 * d.k.y + 3.0 * r2 * d.k.z); + let cross = 2.0 * n.x * n.y * q + 2.0 * d.p.x * n.x + 2.0 * d.p.y * n.y; + return mat2x2( + vec2(radial + 2.0 * n.x * n.x * q + 2.0 * d.p.x * n.y + 6.0 * d.p.y * n.x, cross), + vec2(cross, radial + 2.0 * n.y * n.y * q + 6.0 * d.p.x * n.y + 2.0 * d.p.y * n.x), + ); +} + +// The point distort moves to `goal`, by Newton's method from `goal`. The starting +// error is about k1 r^3, under 1e-3 for usual lenses; each step squares it, so 5 steps +// reach f32's precision. +fn undistort(d: Distortion, goal: vec2) -> vec2 { + var n = goal; + for (var i = 0; i < 5; i++) { + let j = distortion_jacobian(d, n); + let e = distort(d, n) - goal; + let det = j[0].x * j[1].y - j[1].x * j[0].y; + n -= vec2(j[1].y * e.x - j[1].x * e.y, j[0].x * e.y - j[0].y * e.x) / det; + } + return n; +} diff --git a/crates/resolved/src/shaders/sensor.wgsl b/crates/resolved/src/shaders/sensor.wgsl index a4ef28e..7a18fc7 100644 --- a/crates/resolved/src/shaders/sensor.wgsl +++ b/crates/resolved/src/shaders/sensor.wgsl @@ -11,6 +11,8 @@ struct Sensor { detector: vec4, // Width, height, channel count, largest DN. dims: vec4, + // The lens distortion (distortion.wgsl, included first). + distortion: Distortion, } @group(0) @binding(0) var radiance: texture_2d; @@ -28,11 +30,18 @@ fn value(i: u32) -> u32 { let l = textureLoad(radiance, vec2(x, y), 0)[i % n]; // Off-axis falloff: cos^4 of the angle between this pixel's line of sight and the - // boresight. + // boresight. With distortion, the line of sight is the undistorted one, and the + // light spreads over the area the distortion stretches it to. let k = s.intrinsics; - let ray = vec3((f32(x) + 0.5 - k.z) / k.x, (f32(y) + 0.5 - k.w) / k.y, 1.0); + var sight = vec2((f32(x) + 0.5 - k.z) / k.x, (f32(y) + 0.5 - k.w) / k.y); + var stretch = 1.0; + if s.distortion.k.w > 0.5 { + sight = undistort(s.distortion, sight); + stretch = determinant(distortion_jacobian(s.distortion, sight)); + } + let ray = vec3(sight, 1.0); let c2 = 1.0 / dot(ray, ray); - let electrons = clamp(l * s.detector.x * c2 * c2, 0.0, s.detector.y); + let electrons = clamp(l * s.detector.x * c2 * c2 / stretch, 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))); } diff --git a/crates/resolved/tests/gpu/bodies.rs b/crates/resolved/tests/gpu/bodies.rs index 2d0abeb..b52fdb1 100644 --- a/crates/resolved/tests/gpu/bodies.rs +++ b/crates/resolved/tests/gpu/bodies.rs @@ -6,7 +6,7 @@ use crate::common; 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::targets::{FrameTargets, TargetLayout}; use resolved::gpu::GpuContext; use resolved::nalgebra::{UnitQuaternion, Vector3}; use resolved::radiometry::{Radiometry, AU, SUN_RADIUS}; @@ -60,12 +60,7 @@ fn render_limb( ) -> Frame { let targets = FrameTargets::new( ctx, - config.width, - config.height, - config.channels.len() as u32, - 1, - 1, - 0, + &TargetLayout::new(config.width, config.height, config.channels.len() as u32), ); let pass = BodiesPass::new(ctx, 1, false); pass.set_limb_darkening(ctx, limb); diff --git a/crates/resolved/tests/gpu/distortion.rs b/crates/resolved/tests/gpu/distortion.rs new file mode 100644 index 0000000..9cc0005 --- /dev/null +++ b/crates/resolved/tests/gpu/distortion.rs @@ -0,0 +1,340 @@ +//! Camera tier C1's lens distortion: the warp against the projection model, on the GPU +//! and on the CPU (lavapipe). + +use crate::common; +use crate::pixel_area::plate; + +use resolved::gpu::sensor::{SensorPass, SensorUniforms}; +use resolved::gpu::targets::{FrameTargets, TargetLayout}; +use resolved::gpu::AdapterChoice; +use resolved::nalgebra::{UnitQuaternion, Vector2, Vector3}; +use resolved::radiometry::{Radiometry, AU}; +use resolved::{ + BrownConrady, CameraConfig, CameraTier, Fidelity, Frame, Pose, Renderer, RendererOptions, + Scene, Target, +}; + +const ADAPTERS: [AdapterChoice; 2] = [AdapterChoice::Gpu, AdapterChoice::Cpu]; + +/// Strong distortion with all five terms: the corners move about 5 px. +fn lens() -> BrownConrady { + BrownConrady { + k1: 0.3, + k2: -0.2, + k3: 0.1, + p1: 1e-3, + p2: -8e-4, + } +} + +/// The example camera with `distortion`, its principal point off center. +fn camera(distortion: Option) -> CameraConfig { + CameraConfig { + principal_point: Some(Vector2::new(955.25, 543.5)), + distortion, + ..CameraConfig::example() + } +} + +fn renderer(adapter: &AdapterChoice, camera: CameraConfig, tier: CameraTier) -> Renderer { + let options = RendererOptions { + adapter: adapter.clone(), + fidelity: Fidelity { + camera: tier, + ..Default::default() + }, + read_radiance: true, + ..Default::default() + }; + Renderer::new(camera, options).unwrap() +} + +/// A camera 4.2e4 km from Earth on +X, looking along +Z: Earth is 90 degrees off to the +/// side. +fn deep_space() -> Pose { + Pose::new(Vector3::new(4.2e7, 0.0, 0.0), UnitQuaternion::identity()) +} + +/// The Sun behind `camera`, opposite `toward`, lighting what lies that way face on. +fn scene(camera: &Pose, toward: Vector3, moon: Vector3, targets: Vec) -> Scene { + Scene { + earth_attitude: UnitQuaternion::identity(), + sun: camera.position - toward.normalize() * AU, + moon, + moon_attitude: UnitQuaternion::identity(), + targets, + } +} + +fn render(renderer: &mut Renderer, scene: &Scene, camera: &Pose) -> Frame { + renderer.submit(scene, camera).unwrap().wait().unwrap() +} + +/// The green channel of a frame's radiance. +fn green(frame: &Frame) -> Vec { + frame + .radiance + .as_ref() + .unwrap() + .iter() + .map(|p| p[1] as f64) + .collect() +} + +#[test] +fn no_distortion_leaves_the_image_unchanged() { + // Zero coefficients still warp, through a map that moves nothing: the image equals + // the unwarped one. At C0 a lens's distortion is ignored. + let _devices = common::device_lock(); + let view = deep_space(); + for adapter in &ADAPTERS { + let k = camera(None).intrinsics(); + let moon = view.position + k.ray(&Vector2::new(700.3, 400.6)) * 3.844e8; + let mut targets = Vec::new(); + let mut frames = Vec::new(); + for distortion in [None, Some(BrownConrady::default())] { + let mut r = renderer(adapter, camera(distortion), CameraTier::C1); + let model = r.load_model(&plate(1.0)).unwrap(); + targets = vec![Target { + model, + pose: Pose::new( + view.position + Vector3::new(0.031, 0.012, 1.0) * 60.0, + UnitQuaternion::identity(), + ), + }]; + frames.push(render( + &mut r, + &scene(&view, Vector3::z(), moon, targets.clone()), + &view, + )); + } + let (a, b) = (green(&frames[0]), green(&frames[1])); + let peak = a.iter().fold(0.0f64, |m, &x| m.max(x)); + let worst = a + .iter() + .zip(&b) + .fold(0.0f64, |m, (x, y)| m.max((x - y).abs())); + assert!(worst <= 1e-6 * peak, "{adapter:?}: {worst} of {peak}"); + + let mut dn = Vec::new(); + for distortion in [None, Some(lens())] { + let mut r = renderer(adapter, camera(distortion), CameraTier::C0); + r.load_model(&plate(1.0)).unwrap(); + dn.push( + render( + &mut r, + &scene(&view, Vector3::z(), moon, targets.clone()), + &view, + ) + .dn, + ); + } + assert_eq!(dn[0], dn[1], "{adapter:?}"); + } +} + +#[test] +fn points_land_where_the_model_projects_them() { + // A full Moon 4 pixels in radius placed at the center, the edges and the corners, + // with the Sun straight behind the camera from it: the centroid of its light lands + // where Intrinsics::project puts its center, to 0.05 px. At the corners that is over + // 3.5 px from where a pinhole would put it. (With the Sun behind the camera along the + // boresight instead, a Moon 0.1 rad off axis is seen at a phase of 0.1 rad, and its + // light's centroid moves 0.16 px away from the center.) + let _devices = common::device_lock(); + let view = deep_space(); + let k = camera(Some(lens())).intrinsics(); + let pinhole = camera(None).intrinsics(); + for adapter in &ADAPTERS { + let mut r = renderer(adapter, camera(Some(lens())), CameraTier::C1); + let name = r.adapter_info().name.clone(); + for (u, v) in [ + (40.3, 35.7), + (960.2, 30.4), + (1880.6, 40.1), + (25.9, 540.3), + (955.25, 543.5), + (1890.4, 520.8), + (35.2, 1045.5), + (970.7, 1050.2), + (1882.3, 1043.9), + ] { + let uv = Vector2::new(u, v); + let direction = k.ray(&uv); + let image = green(&render( + &mut r, + &scene( + &view, + direction, + view.position + direction * 3.844e9, + Vec::new(), + ), + &view, + )); + let (mut total, mut cu, mut cv) = (0.0, 0.0, 0.0); + for row in (v as usize).saturating_sub(15)..(v as usize + 16).min(1080) { + for col in (u as usize).saturating_sub(15)..(u as usize + 16).min(1920) { + let l = image[row * 1920 + col]; + total += l; + cu += l * (col as f64 + 0.5); + cv += l * (row as f64 + 0.5); + } + } + let centroid = Vector2::new(cu, cv) / total; + assert!( + (centroid - uv).norm() < 0.05, + "{name}: Moon at ({u}, {v}): centroid ({}, {})", + centroid.x, + centroid.y + ); + if !(400.0..1500.0).contains(&u) && !(200.0..900.0).contains(&v) { + let moved = (pinhole.project(&direction).unwrap() - uv).norm(); + assert!(moved > 3.5, "{name}: ({u}, {v}) moved only {moved} px"); + } + } + } +} + +#[test] +fn a_straight_edge_bends_as_the_model_predicts() { + // A plate's straight edge, 0.1 rad right of the boresight, running the height of the + // frame and tilted 0.1 rad from the vertical. In each row its position, the sum of + // the dark fraction of the pixels across it, is compared with where the distorted + // edge crosses that row. A mesh edge's coverage comes in steps, at 4 sample offsets + // across each pixel when rendering 2 x 2 with 2 samples, so a row can be off by up + // to 1/8 px; the tilt sweeps the edge across the offsets, and the error averaged + // over 32 consecutive rows is within 0.05 px. + let _devices = common::device_lock(); + let view = deep_space(); + let k = camera(Some(lens())).intrinsics(); + let d = lens(); + let (a, tilt, range) = (0.1, 0.1f64, 100.0); + // The edge's points in normalized coordinates, along its length t in meters. + let edge = |t: f64| Vector2::new(a - t * tilt.sin() / range, t * tilt.cos() / range); + // Where the distorted edge crosses row v's center line. + let predicted = |v: usize| { + let row = (v as f64 + 0.5 - k.cy) / k.fy; + let mut t = row * range; + for _ in 0..50 { + let h = 1e-6; + let f = d.distort(&edge(t)).y - row; + let slope = (d.distort(&edge(t + h)).y - d.distort(&edge(t - h)).y) / (2.0 * h); + t -= f / slope; + } + k.fx * d.distort(&edge(t)).x + k.cx + }; + let rows: Vec = (16..1064).collect(); + let want: Vec = rows.iter().map(|&v| predicted(v)).collect(); + for adapter in &ADAPTERS { + let mut r = renderer(adapter, camera(Some(lens())), CameraTier::C1); + let name = r.adapter_info().name.clone(); + let model = r.load_model(&plate(10.0)).unwrap(); + let target = Target { + model, + pose: Pose::new( + view.position + Vector3::new(a * range, 0.0, range), + UnitQuaternion::from_axis_angle(&Vector3::z_axis(), tilt), + ), + }; + let moon = Vector3::new(0.0, -3.844e8, 0.0); + let image = green(&render( + &mut r, + &scene(&view, Vector3::z(), moon, vec![target]), + &view, + )); + let errors: Vec = rows + .iter() + .zip(&want) + .map(|(&v, &want)| { + let start = want.floor() as usize - 8; + let inside = image[v * 1920 + start + 24]; + let dark: f64 = (start..start + 16) + .map(|u| 1.0 - image[v * 1920 + u] / inside) + .sum(); + start as f64 + dark - want + }) + .collect(); + let worst_row = errors.iter().fold(0.0f64, |m, e| m.max(e.abs())); + let worst_mean = errors + .windows(32) + .map(|w| (w.iter().sum::() / 32.0).abs()) + .fold(0.0f64, f64::max); + eprintln!( + "straight edge, {name}: worst row off by {worst_row:.4} px, worst 32-row mean by {worst_mean:.4}" + ); + assert!(worst_row < 0.15, "{name}: {worst_row}"); + assert!(worst_mean < 0.05, "{name}: {worst_mean}"); + } +} + +#[test] +fn the_sensor_uses_the_undistorted_line_of_sight() { + // A white surface filling the frame: each pixel's signal falls off with cos^4 of its + // real, undistorted line of sight, and with the area the distortion stretches the + // light over. Against the CPU, to 1 DN. + let config = camera(Some(lens())); + let k = config.intrinsics(); + let d = lens(); + let r = Radiometry::new(&config); + let white: [f32; 4] = std::array::from_fn(|c| { + if c < 3 { + r.lambert_radiance(c, 1.0, 1.0, AU) as f32 + } else { + 0.0 + } + }); + for ctx in common::contexts() { + let targets = FrameTargets::new(ctx, &TargetLayout::new(1920, 1080, 3)); + ctx.queue.write_texture( + wgpu::TexelCopyTextureInfo { + texture: &targets.radiance, + mip_level: 0, + origin: wgpu::Origin3d::ZERO, + aspect: wgpu::TextureAspect::All, + }, + bytemuck::cast_slice(&vec![white; 1920 * 1080]), + wgpu::TexelCopyBufferLayout { + offset: 0, + bytes_per_row: Some(1920 * 16), + rows_per_image: Some(1080), + }, + targets.radiance.size(), + ); + let mut encoder = ctx.device.create_command_encoder(&Default::default()); + SensorPass::new(ctx).encode( + ctx, + &mut encoder, + &targets, + &SensorUniforms::new(&config, &r, true), + None, + ); + ctx.queue.submit([encoder.finish()]); + let dn = targets.read_dn(ctx); + let mut off_by_one = 0; + for v in (0..1080).step_by(37).chain([1079]) { + for u in (0..1920).step_by(37).chain([1919]) { + let uv = Vector2::new(u as f64 + 0.5, v as f64 + 0.5); + let ray = k.ray(&uv); + let n = Vector2::new(ray.x / ray.z, ray.y / ray.z); + let stretch = d.jacobian(&n).determinant(); + let electrons = r.electrons(white[1] as f64, ray.z) / stretch; + let want = config.detector.dn(electrons); + let got = dn[(v * 1920 + u) * 3 + 1]; + assert!( + got.abs_diff(want) <= 1, + "{} ({u}, {v}): {got} vs {want}", + ctx.info.name + ); + off_by_one += usize::from(got != want); + } + } + // The corner reads less than the pinhole's cos^4 alone would give. + let pinhole = camera(None).intrinsics().ray(&Vector2::new(0.5, 0.5)).z; + assert!(dn[1] < config.detector.dn(r.electrons(white[1] as f64, pinhole))); + assert!( + off_by_one < 20, + "{}: {off_by_one} off by one", + ctx.info.name + ); + } +} diff --git a/crates/resolved/tests/gpu/main.rs b/crates/resolved/tests/gpu/main.rs index f20d969..ea7a50e 100644 --- a/crates/resolved/tests/gpu/main.rs +++ b/crates/resolved/tests/gpu/main.rs @@ -6,6 +6,7 @@ mod common; mod analytic; mod bodies; +mod distortion; mod earth_texture; mod earthshine; mod eclipses; diff --git a/crates/resolved/tests/gpu/mesh.rs b/crates/resolved/tests/gpu/mesh.rs index 057e71c..db61bd5 100644 --- a/crates/resolved/tests/gpu/mesh.rs +++ b/crates/resolved/tests/gpu/mesh.rs @@ -10,7 +10,7 @@ use resolved::bodies::{Appearance, EARTH_EQUATORIAL_RADIUS}; 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::targets::{FrameTargets, TargetLayout}; use resolved::gpu::GpuContext; use resolved::model::{LoadOptions, Material, Mesh, ModelId, Part}; use resolved::nalgebra::{UnitQuaternion, Vector3}; @@ -93,12 +93,7 @@ fn render(ctx: &GpuContext, name: &str, scene: &Scene, camera: &Pose, mesh: &Mes let radiometry = Radiometry::new(&config); let targets = FrameTargets::new( ctx, - config.width, - config.height, - config.channels.len() as u32, - 1, - 1, - 0, + &TargetLayout::new(config.width, config.height, config.channels.len() as u32), ); let bodies = BodiesPass::new(ctx, 1, false); let mut meshes = MeshPass::new(ctx, 1); diff --git a/crates/resolved/tests/gpu/pixel_area.rs b/crates/resolved/tests/gpu/pixel_area.rs index d77c6ad..885580b 100644 --- a/crates/resolved/tests/gpu/pixel_area.rs +++ b/crates/resolved/tests/gpu/pixel_area.rs @@ -60,7 +60,7 @@ fn render(renderer: &mut Renderer, scene: &Scene, camera: &Pose) -> Vec { /// A gray square plate in the body XY plane, facing body -Z, from x = 0 to `2 scale` /// and y = `-scale` to `scale`: its left edge is the body Y axis. -fn plate(scale: f64) -> Mesh { +pub fn plate(scale: f64) -> Mesh { let mut s = String::from("solid plate\n"); let c = [[0.0, -1.0], [2.0, -1.0], [2.0, 1.0], [0.0, 1.0]]; for tri in [[c[0], c[1], c[2]], [c[0], c[2], c[3]]] { diff --git a/crates/resolved/tests/gpu/psf.rs b/crates/resolved/tests/gpu/psf.rs index 3bbea1b..c7bd0fa 100644 --- a/crates/resolved/tests/gpu/psf.rs +++ b/crates/resolved/tests/gpu/psf.rs @@ -7,7 +7,7 @@ use crate::common; use resolved::camera::MAX_PSF_SIGMA; use resolved::gpu::blur::{averaged, BlurPass}; -use resolved::gpu::targets::FrameTargets; +use resolved::gpu::targets::{FrameTargets, TargetLayout}; use resolved::gpu::{AdapterChoice, GpuContext}; use resolved::nalgebra::{UnitQuaternion, Vector2, Vector3}; use resolved::radiometry::AU; @@ -35,7 +35,12 @@ fn blur( 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 layout = TargetLayout { + scale, + margin: psf.radius(scale), + ..TargetLayout::new(size.0, size.1, 3) + }; + let targets = FrameTargets::new(ctx, &layout); let (w, h) = (targets.radiance.width(), targets.radiance.height()); assert_eq!(radiance.len(), (w * h) as usize); ctx.queue.write_texture( @@ -53,7 +58,7 @@ fn blur( }, targets.radiance.size(), ); - let pass = BlurPass::new(ctx, psf, 3, scale).unwrap(); + let pass = BlurPass::new(ctx, psf, 3, scale, [size.0, size.1], None).unwrap(); let mut encoder = ctx.device.create_command_encoder(&Default::default()); pass.encode(ctx, &mut encoder, &targets, None); ctx.queue.submit([encoder.finish()]); diff --git a/crates/resolved/tests/gpu/sensor.rs b/crates/resolved/tests/gpu/sensor.rs index 540ee6b..42f17f3 100644 --- a/crates/resolved/tests/gpu/sensor.rs +++ b/crates/resolved/tests/gpu/sensor.rs @@ -6,7 +6,7 @@ use crate::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::targets::{FrameTargets, TargetLayout}; use resolved::gpu::GpuContext; use resolved::nalgebra::{UnitQuaternion, Vector2, Vector3}; use resolved::radiometry::{Radiometry, AU}; @@ -16,7 +16,7 @@ use resolved::{CameraConfig, Pose, Scene, Spectrum}; /// returns the DN image. fn sense(ctx: &GpuContext, name: &str, radiance: &[[f32; 4]]) -> Vec { let config = CameraConfig::example(); - let targets = FrameTargets::new(ctx, config.width, config.height, 3, 1, 1, 0); + let targets = FrameTargets::new(ctx, &TargetLayout::new(config.width, config.height, 3)); ctx.queue.write_texture( wgpu::TexelCopyTextureInfo { texture: &targets.radiance, @@ -33,7 +33,7 @@ fn sense(ctx: &GpuContext, name: &str, radiance: &[[f32; 4]]) -> Vec { targets.radiance.size(), ); let pass = SensorPass::new(ctx); - let uniforms = SensorUniforms::new(&config, &Radiometry::new(&config)); + let uniforms = SensorUniforms::new(&config, &Radiometry::new(&config), false); let mut encoder = ctx.device.create_command_encoder(&Default::default()); pass.encode(ctx, &mut encoder, &targets, &uniforms, None); ctx.queue.submit([encoder.finish()]); @@ -157,9 +157,9 @@ fn sunlit_earth_through_the_whole_chain() { config.near_clip, &body_radiometry, ); - let sensor = SensorUniforms::new(&config, &radiometry); + let sensor = SensorUniforms::new(&config, &radiometry, false); for ctx in common::contexts() { - let targets = FrameTargets::new(ctx, config.width, config.height, 3, 1, 1, 0); + let targets = FrameTargets::new(ctx, &TargetLayout::new(config.width, config.height, 3)); let mut encoder = ctx.device.create_command_encoder(&Default::default()); BodiesPass::new(ctx, 1, false).encode(ctx, &mut encoder, &targets, &bodies, None); SensorPass::new(ctx).encode(ctx, &mut encoder, &targets, &sensor, None); diff --git a/docs/optics.md b/docs/optics.md index 4bdba6f..9f2223f 100644 --- a/docs/optics.md +++ b/docs/optics.md @@ -61,9 +61,82 @@ Some quantities that follow from that: IFOV. Real lenses are not perfect pinholes. They bend straight lines slightly -(distortion), not modeled yet, and spread each point into a small blur (the -point spread function, below). The pinhole model is the starting point -those effects are added to. +(distortion, next) and spread each point into a small blur (the point +spread function, below). The pinhole model is the starting point those +effects are added to. + +### Lens distortion + +Code: `camera.rs` (`BrownConrady`, `Intrinsics`); `gpu/warp.rs`, +`gpu/blur.rs`, `shaders/blur.wgsl`, `shaders/distortion.wgsl`. + +A real lens moves each point a little away from where the pinhole puts it, +by an amount that grows toward the edge of the image. Straight lines in the +scene come out curved. When points move outward more the farther out they +are, the corners stretch (pincushion distortion); when they move inward, +the image bulges (barrel distortion). A lens element slightly off center +adds a smaller sideways shift (tangential distortion). + +The standard model, Brown-Conrady, works on the pinhole's normalized +coordinates x = X / Z and y = Y / Z, with r^2 = x^2 + y^2: + +```text +x_d = x (1 + k1 r^2 + k2 r^4 + k3 r^6) + 2 p1 x y + p2 (r^2 + 2 x^2) +y_d = y (1 + k1 r^2 + k2 r^4 + k3 r^6) + p1 (r^2 + 2 y^2) + 2 p2 x y +u = (f / p) x_d + cx, v = (f / p) y_d + cy +``` + +k1, k2 and k3 are the radial terms, positive for pincushion and negative +for barrel; p1 and p2 are the tangential terms. This is OpenCV's convention, +so the coefficients a camera calibration produces apply directly (OpenCV +lists them in the order k1, k2, p1, p2, k3). For the example camera, k1 = +0.1 alone moves the corners 1.6 px outward. + +That matters for navigation, not just for looks. A limb fit, a star +identification or a target's centroid all assume where each pixel looks; +an unmodeled 1.6 px is a 0.2 mrad pointing error at the corner, which at +36,000 km is 6 km. `Intrinsics::project` and `Intrinsics::ray` include the +distortion (the inverse by Newton's method, to 1e-10 px), so navigation code +and the renderer share one model; the projection matches OpenCV's to 1e-9 +px. + +Rendering it, from camera tier C1 (C0 is always a pinhole): the GPU draws +the target's triangles with a linear projection, which keeps straight edges +straight, and moving their corners alone would leave long edges straight +when they should curve. So the renderer draws the whole frame without +distortion, the bodies and targets sorted together as usual, on a slightly +larger and offset area that covers everything the lens bends into view. +Then it resamples it: for each pixel of the focal plane it finds the +undistorted point the lens bends onto that pixel and reads the render +there, interpolating between the four nearest render pixels. The blur's +first pass (below) does this as it reads its input, so the distorted image +is never written out. Radiance is +the same at both ends of a ray, so moving it is all distortion does to an +extended scene. The resampling works on the finer image C1 renders (below), +before the blur, which is where the lens puts it in a real camera. + +Two effects of distortion fall on the detector instead. The cos^4 falloff +uses each pixel's real line of sight, the undistorted one. And where the +lens stretches the image, the same light covers more pixels: each pixel's +signal is divided by how much the distortion stretches area there (the +determinant of its derivative), 0.6 % at the corners for k1 = 0.1. + +Measured with coefficients of 0.3, -0.2, 0.1 (radial) and 0.001, -0.0008 +(tangential), which move the corners about 4 px: a small Moon's light lands +centered where `Intrinsics::project` puts it to within 0.004 px across the +frame, and a long straight edge follows the predicted curve to within 0.07 +px in every row and 0.02 px averaged over rows. The interpolation adds a +little blur where it falls between render pixels, at most a variance of +1/16 px^2 at 2 x 2 rendering (computed), which takes sigma 0.5 px to 0.56 px +at worst. + +Distortion costs 0.25 to 0.32 ms per 1920 x 1080 frame on an RTX 3060 at +2 x 2 rendering, and 1.1 to 1.3 ms at 3840 x 2160, almost all of it the +blur's first pass reading through it, and under 0.1 ms the larger +undistorted render. A separate pass writing the distorted image out and the +blur reading it back cost 0.7 to 0.8 ms and 2.8 to 3.3 ms: at those sizes +both are limited by how fast the card moves the image through memory. On +lavapipe distortion adds 10 to 11 ms to 1920 x 1080 frames. ### Coordinate conventions @@ -1252,7 +1325,6 @@ velocity as an input, which the scene does not carry yet. ## Not modeled yet - Earthshine on the Moon's night side -- Lens distortion - The diffraction pattern's rings, and a lens's measured point spread function - Noise: photon shot noise, read noise, dark current, pixel-to-pixel variation - Saturation spilling into neighbors (blooming) diff --git a/readme.md b/readme.md index 021f39b..8ca358f 100644 --- a/readme.md +++ b/readme.md @@ -189,10 +189,11 @@ and targets to DN through `Renderer`, and the CLI renders scene files. Lighting L1 draws Earth's surface from Blue Marble and the Moon's from the LRO Camera's map, shadows on the target, light from Earth and the Moon on it, eclipses on Earth and the Moon, and the Sun's limb darkening. Camera C1 is under way: so -far, pixels that collect light over their area, and lens blur (a Gaussian -PSF) applied to the scene rendered 2 x 2 finer than the pixels. Planned, in order: +far, pixels that collect light over their area, lens distortion +(Brown-Conrady), and lens blur (a Gaussian PSF) applied to the scene +rendered 2 x 2 finer than the pixels. Planned, in order: -- Camera: lens distortion, apparent positions, stars, +- Camera: apparent positions, stars, detector noise, pixel-to-pixel variation, saturation and blooming. - Lighting: BRDF and spectral materials, atmospheric scattering, clouds, lunar terrain.