Added spectrum and camera specs

crates/resolved-cli/Cargo.toml

@@ -12,6 +12,8 @@ repository.workspace = true
 [[bin]]
 name = "resolved"
 path = "src/main.rs"
+# The library is also named `resolved`; documenting both collides in target/doc.
+doc = false
 
 [dependencies]
 resolved.workspace = true

crates/resolved/Cargo.toml

@@ -13,6 +13,7 @@ repository.workspace = true
 
 [dependencies]
 nalgebra = "0.35"
+thiserror = "2.0"
 
 [dev-dependencies]
 approx = "0.5"

crates/resolved/src/camera.rs

@@ -0,0 +1,281 @@
+//! Camera geometry, optics and detector parameters.
+//!
+//! Pixel coordinates are continuous: pixel `(i, j)` covers `[i, i + 1) x [j, j + 1)`, so
+//! its center is at `(i + 0.5, j + 0.5)`. `u` increases with column and `v` with row,
+//! 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 crate::scene::Pose;
+use crate::spectrum::Spectrum;
+use crate::Error;
+
+/// One color channel of the detector, such as the red pixels behind a Bayer filter.
+#[derive(Clone, Debug, PartialEq)]
+pub struct Channel {
+    /// A label for output files and logs.
+    pub name: String,
+    /// Electrons per incident photon, including any filter in front of these pixels.
+    pub quantum_efficiency: Spectrum,
+}
+
+/// How the detector turns electrons into numbers.
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub struct Detector {
+    /// Electrons per DN.
+    pub gain: f64,
+    /// DN added to every pixel before quantization.
+    pub bias: f64,
+    /// Electrons a pixel holds before it saturates.
+    pub full_well: f64,
+    /// Bits per pixel out of the ADC; output saturates at `2^adc_bits - 1`.
+    pub adc_bits: u32,
+}
+
+/// A pinhole camera with radiometric optics and a detector.
+#[derive(Clone, Debug, PartialEq)]
+pub struct CameraConfig {
+    /// Image width in pixels.
+    pub width: u32,
+    /// Image height in pixels.
+    pub height: u32,
+    /// Distance from the pinhole to the focal plane, in meters.
+    pub focal_length: f64,
+    /// Center-to-center pixel spacing, in meters. Pixels are square.
+    pub pixel_pitch: f64,
+    /// Where the boresight meets the image, in pixels. `None` puts it at the center.
+    pub principal_point: Option<Vector2<f64>>,
+    /// Focal length over entrance pupil diameter.
+    pub f_number: f64,
+    /// Integration time, in seconds.
+    pub exposure: f64,
+    /// Fraction of light the optics pass, excluding per-channel filters.
+    pub transmission: Spectrum,
+    /// The detector's color channels.
+    pub channels: Vec<Channel>,
+    /// Electrons to DN.
+    pub detector: Detector,
+    /// Nearest distance that is drawn, in meters.
+    pub near_clip: f64,
+    /// The camera frame in the observer's body frame.
+    pub mount: Pose,
+}
+
+impl CameraConfig {
+    /// A documented example camera, small enough to check its numbers by hand.
+    ///
+    /// - 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.
+    /// - 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.
+    /// - 0.1 ms exposure: a white Lambertian surface facing the Sun fills about half the
+    ///   well (44-54 % across the channels), so Earth and the target are both unsaturated.
+    /// - Mounted with no offset or rotation on the observer.
+    pub fn example() -> Self {
+        let channel = |name: &str, center: f64| Channel {
+            name: name.to_owned(),
+            quantum_efficiency: Spectrum::gaussian(0.6, center, 80e-9),
+        };
+        Self {
+            width: 1920,
+            height: 1080,
+            focal_length: 50e-3,
+            pixel_pitch: 5.5e-6,
+            principal_point: None,
+            f_number: 4.0,
+            exposure: 1e-4,
+            transmission: Spectrum::constant(0.9),
+            channels: vec![
+                channel("r", 620e-9),
+                channel("g", 540e-9),
+                channel("b", 460e-9),
+            ],
+            detector: Detector {
+                gain: 5.0,
+                bias: 64.0,
+                full_well: 20_000.0,
+                adc_bits: 12,
+            },
+            near_clip: 0.01,
+            mount: Pose::new(Vector3::zeros(), UnitQuaternion::identity()),
+        }
+    }
+
+    /// Checks that every parameter describes something physically possible.
+    pub fn validate(&self) -> Result<(), Error> {
+        let fail = |what: &str| Err(Error::InvalidCamera(what.to_owned()));
+        let positive = |x: f64| x.is_finite() && x > 0.0;
+        if self.width == 0 || self.height == 0 {
+            return fail("image size must be nonzero");
+        }
+        if !positive(self.focal_length) || !positive(self.pixel_pitch) {
+            return fail("focal length and pixel pitch must be positive");
+        }
+        if !positive(self.f_number) || !positive(self.exposure) || !positive(self.near_clip) {
+            return fail("f-number, exposure and near clip must be positive");
+        }
+        if let Some(p) = self.principal_point {
+            if !p.iter().all(|x| x.is_finite()) {
+                return fail("principal point must be finite");
+            }
+        }
+        if self.channels.is_empty() {
+            return fail("at least one channel is required");
+        }
+        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");
+        }
+        if !(1..=16).contains(&d.adc_bits) {
+            return fail("ADC bits must be between 1 and 16");
+        }
+        Ok(())
+    }
+
+    /// Pinhole projection parameters in pixels.
+    pub fn intrinsics(&self) -> Intrinsics {
+        let f = self.focal_length / self.pixel_pitch;
+        let c = self
+            .principal_point
+            .unwrap_or_else(|| Vector2::new(self.width as f64, self.height as f64) / 2.0);
+        Intrinsics {
+            fx: f,
+            fy: f,
+            cx: c.x,
+            cy: c.y,
+            width: self.width,
+            height: self.height,
+        }
+    }
+}
+
+/// Pinhole projection in pixels: `u = fx x / z + cx`, `v = fy y / z + cy`.
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub struct Intrinsics {
+    /// Focal length in pixels along `u`.
+    pub fx: f64,
+    /// Focal length in pixels along `v`.
+    pub fy: f64,
+    /// Principal point column, in pixels.
+    pub cx: f64,
+    /// Principal point row, in pixels.
+    pub cy: f64,
+    /// Image width in pixels.
+    pub width: u32,
+    /// Image height in pixels.
+    pub height: u32,
+}
+
+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<f64>) -> Option<Vector2<f64>> {
+        (p.z > 0.0)
+            .then(|| Vector2::new(self.fx * p.x / p.z + self.cx, self.fy * p.y / p.z + self.cy))
+    }
+
+    /// Unit direction in the camera frame of the ray through pixel coordinates `uv`.
+    pub fn ray(&self, uv: &Vector2<f64>) -> Vector3<f64> {
+        Vector3::new((uv.x - self.cx) / self.fx, (uv.y - self.cy) / self.fy, 1.0).normalize()
+    }
+
+    /// Camera frame to wgpu clip space, with reversed-Z depth and no far plane.
+    ///
+    /// 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
+    /// `y` map pixel coordinates onto normalized device coordinates, with row 0 at the
+    /// top (`y = +1`).
+    pub fn clip_from_camera(&self, near: f64) -> Matrix4<f64> {
+        let (w, h) = (self.width as f64, self.height as f64);
+        #[rustfmt::skip]
+        let m = Matrix4::new(
+            2.0 * self.fx / w, 0.0,               2.0 * self.cx / w - 1.0, 0.0,
+            0.0,               -2.0 * self.fy / h, 1.0 - 2.0 * self.cy / h, 0.0,
+            0.0,               0.0,               0.0,                     near,
+            0.0,               0.0,               1.0,                     0.0,
+        );
+        m
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use approx::assert_relative_eq;
+    use nalgebra::Vector4;
+
+    #[test]
+    fn example_camera_is_valid() {
+        let camera = CameraConfig::example();
+        camera.validate().unwrap();
+        let k = camera.intrinsics();
+        let hfov = 2.0 * (k.width as f64 / 2.0 / k.fx).atan();
+        assert_relative_eq!(hfov.to_degrees(), 12.056, epsilon = 1e-3);
+        assert_relative_eq!(1.0 / k.fx, 110e-6, epsilon = 1e-12);
+    }
+
+    #[test]
+    fn validate_rejects_nonphysical_values() {
+        let mut camera = CameraConfig::example();
+        camera.f_number = 0.0;
+        assert!(camera.validate().is_err());
+
+        let mut camera = CameraConfig::example();
+        camera.channels.clear();
+        assert!(camera.validate().is_err());
+
+        let mut camera = CameraConfig::example();
+        camera.detector.adc_bits = 17;
+        assert!(camera.validate().is_err());
+    }
+
+    #[test]
+    fn boresight_projects_to_principal_point() {
+        let k = CameraConfig::example().intrinsics();
+        let uv = k.project(&Vector3::new(0.0, 0.0, 10.0)).unwrap();
+        assert_relative_eq!(uv, Vector2::new(960.0, 540.0), epsilon = 1e-12);
+        assert!(k.project(&Vector3::new(0.0, 0.0, -10.0)).is_none());
+    }
+
+    #[test]
+    fn ray_inverts_project() {
+        let mut camera = CameraConfig::example();
+        camera.principal_point = Some(Vector2::new(955.25, 543.5));
+        let k = camera.intrinsics();
+        for uv in [
+            Vector2::new(0.5, 0.5),
+            Vector2::new(1919.5, 7.25),
+            Vector2::new(400.0, 1000.0),
+        ] {
+            let p = k.ray(&uv) * 37.0;
+            assert_relative_eq!(k.project(&p).unwrap(), uv, epsilon = 1e-9);
+        }
+    }
+
+    #[test]
+    fn clip_matrix_agrees_with_projection() {
+        let mut camera = CameraConfig::example();
+        camera.principal_point = Some(Vector2::new(955.25, 543.5));
+        let k = camera.intrinsics();
+        let m = k.clip_from_camera(camera.near_clip);
+        for p in [
+            Vector3::new(0.3, -0.2, 10.0),
+            Vector3::new(-1.0e5, 4.0e5, 7.0e6),
+            Vector3::new(2.0, 1.0, 1.5e11),
+        ] {
+            let clip = m * Vector4::new(p.x, p.y, p.z, 1.0);
+            let ndc = clip.xyz() / clip.w;
+            let uv = k.project(&p).unwrap();
+            let (w, h) = (k.width as f64, k.height as f64);
+            assert_relative_eq!(ndc.x, 2.0 * uv.x / w - 1.0, epsilon = 1e-12);
+            assert_relative_eq!(ndc.y, 1.0 - 2.0 * uv.y / h, epsilon = 1e-12);
+            assert_relative_eq!(ndc.z, camera.near_clip / p.z, max_relative = 1e-12);
+        }
+        // The near plane sits at depth 1 and depth falls with distance.
+        let near = m * Vector4::new(0.0, 0.0, camera.near_clip, 1.0);
+        assert_relative_eq!(near.z / near.w, 1.0, epsilon = 1e-12);
+    }
+}

crates/resolved/src/error.rs

@@ -0,0 +1,10 @@
+//! The crate's error type.
+
+/// Everything that can go wrong in this crate.
+#[derive(Debug, thiserror::Error)]
+#[non_exhaustive]
+pub enum Error {
+    /// A camera configuration that cannot describe a real camera.
+    #[error("invalid camera configuration: {0}")]
+    InvalidCamera(String),
+}

crates/resolved/src/lib.rs

@@ -13,6 +13,12 @@
 /// version of `Vector3` and `UnitQuaternion` that the API was compiled against.
 pub use nalgebra;
 
+pub mod camera;
+pub mod error;
 pub mod scene;
+pub mod spectrum;
 
+pub use camera::{CameraConfig, Channel, Detector, Intrinsics};
+pub use error::Error;
 pub use scene::{CameraRelative, Placement, Pose, Scene, Target};
+pub use spectrum::Spectrum;

crates/resolved/src/spectrum.rs

@@ -0,0 +1,112 @@
+//! Quantities that vary with wavelength: QE, transmission, reflectance, irradiance.
+//!
+//! Wavelengths are in meters, like every other length in this crate.
+
+/// A function of wavelength, linearly interpolated between samples and zero outside them.
+#[derive(Clone, Debug, PartialEq)]
+pub struct Spectrum(Repr);
+
+#[derive(Clone, Debug, PartialEq)]
+enum Repr {
+    Constant(f64),
+    Sampled {
+        wavelength: Vec<f64>,
+        value: Vec<f64>,
+    },
+}
+
+impl Spectrum {
+    /// The same value at every wavelength.
+    pub fn constant(value: f64) -> Self {
+        Self(Repr::Constant(value))
+    }
+
+    /// Samples at strictly increasing wavelengths, in meters.
+    ///
+    /// Returns `None` if the lengths differ, there are fewer than two samples, or the
+    /// wavelengths are not strictly increasing and finite.
+    pub fn sampled(wavelength: Vec<f64>, value: Vec<f64>) -> Option<Self> {
+        let increasing = wavelength.windows(2).all(|w| w[0] < w[1]);
+        let finite = wavelength.iter().chain(&value).all(|x| x.is_finite());
+        (wavelength.len() == value.len() && wavelength.len() >= 2 && increasing && finite)
+            .then_some(Self(Repr::Sampled { wavelength, value }))
+    }
+
+    /// A Gaussian band with the given peak, center and full width at half maximum, sampled
+    /// out to three widths either side.
+    pub fn gaussian(peak: f64, center: f64, fwhm: f64) -> Self {
+        const SAMPLES: usize = 121;
+        let sigma = fwhm / (8.0 * std::f64::consts::LN_2).sqrt();
+        let (lo, hi) = (center - 3.0 * fwhm, center + 3.0 * fwhm);
+        let wavelength: Vec<f64> = (0..SAMPLES)
+            .map(|i| lo + (hi - lo) * i as f64 / (SAMPLES - 1) as f64)
+            .collect();
+        let value = wavelength
+            .iter()
+            .map(|l| peak * (-0.5 * ((l - center) / sigma).powi(2)).exp())
+            .collect();
+        Self(Repr::Sampled { wavelength, value })
+    }
+
+    /// The value at `wavelength`, in meters.
+    pub fn at(&self, wavelength: f64) -> f64 {
+        match &self.0 {
+            Repr::Constant(v) => *v,
+            Repr::Sampled {
+                wavelength: w,
+                value: v,
+            } => {
+                if wavelength < w[0] || wavelength > w[w.len() - 1] {
+                    return 0.0;
+                }
+                let i = w
+                    .partition_point(|&x| x <= wavelength)
+                    .clamp(1, w.len() - 1);
+                let t = (wavelength - w[i - 1]) / (w[i] - w[i - 1]);
+                v[i - 1] + t * (v[i] - v[i - 1])
+            }
+        }
+    }
+
+    /// The wavelength range outside which this spectrum is zero, or `None` if it is
+    /// constant and so has no bound of its own.
+    pub fn support(&self) -> Option<(f64, f64)> {
+        match &self.0 {
+            Repr::Constant(_) => None,
+            Repr::Sampled { wavelength: w, .. } => Some((w[0], w[w.len() - 1])),
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use approx::assert_relative_eq;
+
+    #[test]
+    fn sampled_interpolates_linearly_and_is_zero_outside() {
+        let s = Spectrum::sampled(vec![400e-9, 500e-9, 600e-9], vec![0.0, 1.0, 0.5]).unwrap();
+        assert_relative_eq!(s.at(450e-9), 0.5, epsilon = 1e-12);
+        assert_relative_eq!(s.at(500e-9), 1.0, epsilon = 1e-12);
+        assert_relative_eq!(s.at(550e-9), 0.75, epsilon = 1e-12);
+        assert_relative_eq!(s.at(600e-9), 0.5, epsilon = 1e-12);
+        assert_eq!(s.at(399e-9), 0.0);
+        assert_eq!(s.at(601e-9), 0.0);
+    }
+
+    #[test]
+    fn sampled_rejects_bad_input() {
+        assert!(Spectrum::sampled(vec![1.0], vec![1.0]).is_none());
+        assert!(Spectrum::sampled(vec![1.0, 2.0], vec![1.0]).is_none());
+        assert!(Spectrum::sampled(vec![2.0, 1.0], vec![1.0, 1.0]).is_none());
+        assert!(Spectrum::sampled(vec![1.0, f64::NAN], vec![1.0, 1.0]).is_none());
+    }
+
+    #[test]
+    fn gaussian_half_maximum_at_half_width() {
+        let s = Spectrum::gaussian(0.6, 540e-9, 80e-9);
+        assert_relative_eq!(s.at(540e-9), 0.6, epsilon = 1e-12);
+        assert_relative_eq!(s.at(500e-9), 0.3, max_relative = 1e-3);
+        assert_relative_eq!(s.at(580e-9), 0.3, max_relative = 1e-3);
+    }
+}