368 lines
14 KiB
Rust
368 lines
14 KiB
Rust
//! Camera tier C1's detector noise against its statistics: the photon transfer curve,
|
|
//! dark frames, exact Poisson counts, reproducibility, and a faint star's SNR, on the
|
|
//! GPU and on the CPU (lavapipe).
|
|
|
|
use crate::common;
|
|
|
|
use std::sync::Arc;
|
|
|
|
use resolved::camera::Detector;
|
|
use resolved::gpu::sensor::{SensorPass, SensorPattern, SensorUniforms};
|
|
use resolved::gpu::targets::{FrameTargets, TargetLayout};
|
|
use resolved::gpu::{AdapterChoice, GpuContext};
|
|
use resolved::nalgebra::{UnitQuaternion, Vector2, Vector3};
|
|
use resolved::radiometry::{Radiometry, AU};
|
|
use resolved::stars::CatalogStar;
|
|
use resolved::{
|
|
CameraConfig, CameraTier, EarthTexture, Fidelity, LightingTier, Pose, Renderer,
|
|
RendererOptions, Scene, StarCatalog, StarPhotometry,
|
|
};
|
|
|
|
const SIZE: u32 = 256;
|
|
|
|
/// 256 x 256 pixels of a 1 m focal length, so the falloff off axis is under 1e-5, with
|
|
/// `detector`.
|
|
fn camera(detector: Detector) -> CameraConfig {
|
|
CameraConfig {
|
|
width: SIZE,
|
|
height: SIZE,
|
|
focal_length: 1.0,
|
|
detector,
|
|
..CameraConfig::example()
|
|
}
|
|
}
|
|
|
|
/// The example detector with its noise replaced and no fixed pattern.
|
|
fn detector(read_noise: f64, dark_current: f64, shot_noise: bool) -> Detector {
|
|
Detector {
|
|
read_noise,
|
|
dark_current,
|
|
shot_noise,
|
|
..CameraConfig::example().detector.noiseless()
|
|
}
|
|
}
|
|
|
|
/// A detector that reads whole electrons: 1 e-/DN, 16 bits, a bias of 100 DN.
|
|
fn counting(read_noise: f64, dark_current: f64) -> Detector {
|
|
Detector {
|
|
gain: 1.0,
|
|
bias: 100.0,
|
|
full_well: 60_000.0,
|
|
adc_bits: 16,
|
|
read_noise,
|
|
dark_current,
|
|
shot_noise: true,
|
|
..CameraConfig::example().detector.noiseless()
|
|
}
|
|
}
|
|
|
|
/// Runs the sensor pass with noise over a flat field that gives `electrons` in every
|
|
/// pixel and channel, for frame `frame` of seed `seed`, and returns the DN image.
|
|
fn flat(
|
|
ctx: &GpuContext,
|
|
config: &CameraConfig,
|
|
electrons: f64,
|
|
seed: u64,
|
|
frame: u64,
|
|
) -> Vec<u16> {
|
|
let targets = FrameTargets::new(ctx, &TargetLayout::new(SIZE, SIZE, 3));
|
|
let radiometry = Radiometry::new(config);
|
|
let level = (electrons / radiometry.electrons_per_radiance) as f32;
|
|
let radiance = vec![[level, level, level, 0.0]; (SIZE * SIZE) as usize];
|
|
ctx.queue.write_texture(
|
|
wgpu::TexelCopyTextureInfo {
|
|
texture: &targets.radiance,
|
|
mip_level: 0,
|
|
origin: wgpu::Origin3d::ZERO,
|
|
aspect: wgpu::TextureAspect::All,
|
|
},
|
|
bytemuck::cast_slice(&radiance),
|
|
wgpu::TexelCopyBufferLayout {
|
|
offset: 0,
|
|
bytes_per_row: Some(SIZE * 16),
|
|
rows_per_image: Some(SIZE),
|
|
},
|
|
targets.radiance.size(),
|
|
);
|
|
let pass = SensorPass::new(ctx, false, true, SensorPattern::None);
|
|
let mut uniforms = SensorUniforms::new(config, &radiometry, false);
|
|
uniforms.set_key(seed, frame);
|
|
let mut encoder = ctx.device.create_command_encoder(&Default::default());
|
|
pass.encode(ctx, &mut encoder, &targets, &uniforms, None);
|
|
ctx.queue.submit([encoder.finish()]);
|
|
targets.read_dn(ctx)
|
|
}
|
|
|
|
fn mean(x: &[f64]) -> f64 {
|
|
x.iter().sum::<f64>() / x.len() as f64
|
|
}
|
|
|
|
fn variance(x: &[f64]) -> f64 {
|
|
let m = mean(x);
|
|
x.iter().map(|v| (v - m) * (v - m)).sum::<f64>() / (x.len() - 1) as f64
|
|
}
|
|
|
|
fn correlation(a: &[f64], b: &[f64]) -> f64 {
|
|
let (ma, mb) = (mean(a), mean(b));
|
|
let cov: f64 = a.iter().zip(b).map(|(x, y)| (x - ma) * (y - mb)).sum();
|
|
cov / (variance(a) * variance(b)).sqrt() / (a.len() - 1) as f64
|
|
}
|
|
|
|
/// Channel `k`'s values of a DN image, as f64.
|
|
fn channel(dn: &[u16], k: usize) -> Vec<f64> {
|
|
dn.iter().skip(k).step_by(3).map(|&d| d as f64).collect()
|
|
}
|
|
|
|
#[test]
|
|
fn flat_fields_follow_the_photon_transfer_curve() {
|
|
// Two frames of a flat field at each of six levels, as EMVA 1288 measures a camera:
|
|
// the mean of the two less the bias, and half the variance of their difference,
|
|
// which leaves out anything the same in both. The temporal variance in DN^2 is
|
|
// (electrons + read noise^2) / gain^2 + 1/12 for the rounding, and the GPU and
|
|
// lavapipe each match it to 3 % (65,536 pixels give a sampling error of 0.6 %),
|
|
// and a line through the levels has slope 1 / gain.
|
|
let config = camera(detector(5.0, 0.0, true));
|
|
let d = config.detector;
|
|
for ctx in common::contexts() {
|
|
let name = &ctx.info.name;
|
|
let mut points = Vec::new();
|
|
for (level, electrons) in [100.0, 500.0, 2000.0, 5000.0, 10_000.0, 16_000.0]
|
|
.into_iter()
|
|
.enumerate()
|
|
{
|
|
// Frames of their own at each level, so the levels sample independently.
|
|
let a = flat(ctx, &config, electrons, 7, 2 * level as u64);
|
|
let b = flat(ctx, &config, electrons, 7, 2 * level as u64 + 1);
|
|
for k in 0..3 {
|
|
let (a, b) = (channel(&a, k), channel(&b, k));
|
|
let m = mean(
|
|
&a.iter()
|
|
.zip(&b)
|
|
.map(|(x, y)| (x + y) / 2.0)
|
|
.collect::<Vec<_>>(),
|
|
) - d.bias;
|
|
let v = variance(&a.iter().zip(&b).map(|(x, y)| x - y).collect::<Vec<_>>()) / 2.0;
|
|
let want = d.variance(electrons, config.exposure) / (d.gain * d.gain);
|
|
assert!(
|
|
(m / (electrons / d.gain) - 1.0).abs() < 2e-3,
|
|
"{name}: {electrons} e-, channel {k}: mean {m} DN"
|
|
);
|
|
assert!(
|
|
(v / want - 1.0).abs() < 0.03,
|
|
"{name}: {electrons} e-, channel {k}: variance {v} DN^2, expected {want}"
|
|
);
|
|
points.push((m, v));
|
|
}
|
|
}
|
|
let (mx, my) = (
|
|
mean(&points.iter().map(|p| p.0).collect::<Vec<_>>()),
|
|
mean(&points.iter().map(|p| p.1).collect::<Vec<_>>()),
|
|
);
|
|
let slope = points.iter().map(|p| (p.0 - mx) * (p.1 - my)).sum::<f64>()
|
|
/ points.iter().map(|p| (p.0 - mx) * (p.0 - mx)).sum::<f64>();
|
|
assert!(
|
|
(slope * d.gain - 1.0).abs() < 0.02,
|
|
"{name}: slope {slope} DN, expected {}",
|
|
1.0 / d.gain
|
|
);
|
|
}
|
|
}
|
|
|
|
#[test]
|
|
fn dark_frames_carry_dark_current_and_read_noise() {
|
|
// No light, 1000 e-/s of dark current over 0.1 s: a mean of 100 electrons, 20 DN
|
|
// above the bias, and a variance of (100 + 5^2) / 25 + 1/12 DN^2. Without shot
|
|
// noise the dark electrons are only an offset, and the variance is the read
|
|
// noise's and the rounding's.
|
|
for ctx in common::contexts() {
|
|
let name = &ctx.info.name;
|
|
for shot in [true, false] {
|
|
let config = CameraConfig {
|
|
exposure: 0.1,
|
|
..camera(detector(5.0, 1000.0, shot))
|
|
};
|
|
let d = config.detector;
|
|
let dn = channel(&flat(ctx, &config, 0.0, 3, 0), 1);
|
|
let m = mean(&dn) - d.bias;
|
|
let v = variance(&dn);
|
|
let want = d.variance(0.0, config.exposure) / (d.gain * d.gain);
|
|
assert!((m - 20.0).abs() < 0.05, "{name}: shot {shot}: mean {m} DN");
|
|
assert!(
|
|
(v / want - 1.0).abs() < 0.03,
|
|
"{name}: shot {shot}: variance {v} DN^2, expected {want}"
|
|
);
|
|
}
|
|
}
|
|
}
|
|
|
|
#[test]
|
|
fn small_counts_are_poisson() {
|
|
// At 1 e-/DN with no read noise the DN less the bias is the count itself. Below a
|
|
// mean of 50 it is drawn exactly: every count's frequency over 196,608 values
|
|
// matches the Poisson probability to 5 standard errors. Above 50 the Gaussian of
|
|
// the same variance takes over: its mean and variance still match.
|
|
for ctx in common::contexts() {
|
|
let name = &ctx.info.name;
|
|
let config = camera(counting(0.0, 0.0));
|
|
for mean_count in [0.5, 3.0, 20.0, 45.0, 60.0] {
|
|
let dn = flat(ctx, &config, mean_count, 11, 0);
|
|
let counts: Vec<f64> = dn.iter().map(|&d| d as f64 - 100.0).collect();
|
|
let n = counts.len() as f64;
|
|
let (m, v) = (mean(&counts), variance(&counts));
|
|
assert!(
|
|
(m / mean_count - 1.0).abs() < 0.02 && (v / mean_count - 1.0).abs() < 0.03,
|
|
"{name}: mean {mean_count}: measured mean {m}, variance {v}"
|
|
);
|
|
if mean_count > 50.0 {
|
|
continue;
|
|
}
|
|
let mut p = (-mean_count).exp();
|
|
for k in 0..(mean_count * 3.0 + 10.0) as usize {
|
|
if k > 0 {
|
|
p *= mean_count / k as f64;
|
|
}
|
|
let seen = counts.iter().filter(|&&c| c == k as f64).count() as f64 / n;
|
|
let error = (p * (1.0 - p) / n).sqrt();
|
|
assert!(
|
|
(seen - p).abs() < 5.0 * error + 1e-5,
|
|
"{name}: mean {mean_count}: P({k}) {seen}, Poisson {p}"
|
|
);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
#[test]
|
|
fn noise_repeats_by_seed_and_frame_and_is_independent() {
|
|
// The same seed and frame give the same image; another frame or another seed gives
|
|
// noise uncorrelated with it (|r| < 0.01 over 196,608 values, against a sampling
|
|
// error of 0.0023), and neighboring values within a frame are uncorrelated too.
|
|
let config = camera(detector(5.0, 0.0, true));
|
|
for ctx in common::contexts() {
|
|
let name = &ctx.info.name;
|
|
let base = flat(ctx, &config, 1000.0, 1, 0);
|
|
assert_eq!(base, flat(ctx, &config, 1000.0, 1, 0), "{name}: repeat");
|
|
let values: Vec<f64> = base.iter().map(|&d| d as f64).collect();
|
|
for (what, other) in [
|
|
("next frame", flat(ctx, &config, 1000.0, 1, 1)),
|
|
("next seed", flat(ctx, &config, 1000.0, 2, 0)),
|
|
(
|
|
"seed's high word",
|
|
flat(ctx, &config, 1000.0, 1 + (1 << 32), 0),
|
|
),
|
|
] {
|
|
let other: Vec<f64> = other.iter().map(|&d| d as f64).collect();
|
|
let r = correlation(&values, &other);
|
|
assert!(r.abs() < 0.01, "{name}: {what}: r = {r}");
|
|
}
|
|
let r = correlation(&values[..values.len() - 1], &values[1..]);
|
|
assert!(r.abs() < 0.01, "{name}: neighboring values: r = {r}");
|
|
let r = correlation(&values[..values.len() - 3], &values[3..]);
|
|
assert!(r.abs() < 0.01, "{name}: neighboring pixels: r = {r}");
|
|
}
|
|
}
|
|
|
|
#[test]
|
|
fn a_faint_star_has_the_signal_to_noise_the_detector_gives() {
|
|
// 400 identical stars of V 11.5 in a 20 x 20 grid, 12 pixels apart, through the
|
|
// example detector's noise, exposed to give each 300 electrons. Each star's electrons summed over the 7 x 7
|
|
// pixels around it, less the bias, average S + 49 I_d t, S its electrons and I_d t
|
|
// each pixel's dark electrons, and vary from star to star with variance
|
|
// S + 49 (I_d t + read^2 + gain^2 / 12) electrons^2: an SNR of about 7. The sums'
|
|
// mean and spread match that to 2 % and 12 % (400 stars measure a spread to 3.5 %). The
|
|
// renderer numbers its frames, and setting the number back repeats a frame.
|
|
let _devices = common::device_lock();
|
|
// The example detector's noise, without its fixed pattern.
|
|
let noisy = detector(5.0, 20.0, true);
|
|
let photometry = StarPhotometry::new(&camera(noisy));
|
|
let signal = 300.0;
|
|
let exposure = signal / (photometry.rate(1, 11.5, 0.6) * photometry.aperture);
|
|
let config = CameraConfig {
|
|
exposure,
|
|
..camera(noisy)
|
|
};
|
|
let k = config.intrinsics();
|
|
let position = Vector3::new(4.2e7, 0.0, 0.0);
|
|
let look = Pose::looking_at(position, position + Vector3::y(), Vector3::z()).unwrap();
|
|
let mut stars = Vec::new();
|
|
let mut centers = Vec::new();
|
|
for j in 0..20 {
|
|
for i in 0..20 {
|
|
let uv = Vector2::new(14.3 + 12.0 * i as f64, 14.7 + 12.0 * j as f64);
|
|
let d = look.attitude * k.ray(&uv).normalize();
|
|
stars.push(CatalogStar {
|
|
ra: d.y.atan2(d.x),
|
|
dec: d.z.asin(),
|
|
pm_ra: 0.0,
|
|
pm_dec: 0.0,
|
|
v: 11.5,
|
|
bv: 0.6,
|
|
});
|
|
centers.push(uv);
|
|
}
|
|
}
|
|
let d = config.detector;
|
|
let dark = d.dark_current * exposure;
|
|
let per_pixel = d.variance(0.0, exposure);
|
|
let spread = (signal + 49.0 * per_pixel).sqrt();
|
|
let scene = Scene {
|
|
earth_attitude: UnitQuaternion::identity(),
|
|
sun: position - Vector3::y() * AU,
|
|
moon: Vector3::new(0.0, 0.0, -3.844e8),
|
|
moon_attitude: UnitQuaternion::identity(),
|
|
targets: Vec::new(),
|
|
velocities: Default::default(),
|
|
};
|
|
let catalog = Arc::new(StarCatalog::new(stars, 0.0));
|
|
for adapter in [AdapterChoice::Gpu, AdapterChoice::Cpu] {
|
|
let black = EarthTexture::from_rgba8(1, 1, vec![0, 0, 0, 255]).unwrap();
|
|
let options = RendererOptions {
|
|
adapter,
|
|
fidelity: Fidelity {
|
|
lighting: LightingTier::L1,
|
|
camera: CameraTier::C1,
|
|
},
|
|
earth_texture: Some(Arc::new(black)),
|
|
moon_texture: Some(common::uniform_moon()),
|
|
stars: Some(catalog.clone()),
|
|
noise_seed: 5,
|
|
..Default::default()
|
|
};
|
|
let mut r = Renderer::new(config.clone(), options).unwrap();
|
|
let name = r.adapter_info().name.clone();
|
|
let first = r.submit(&scene, &look).unwrap().wait().unwrap();
|
|
let second = r.submit(&scene, &look).unwrap().wait().unwrap();
|
|
assert_eq!((first.number, second.number), (0, 1), "{name}");
|
|
assert_ne!(first.dn, second.dn, "{name}");
|
|
r.set_frame_number(0);
|
|
let again = r.submit(&scene, &look).unwrap().wait().unwrap();
|
|
assert_eq!(again.number, 0, "{name}");
|
|
assert_eq!(first.dn, again.dn, "{name}: frame 0 again");
|
|
|
|
let sums: Vec<f64> = centers
|
|
.iter()
|
|
.map(|c| {
|
|
let (u0, v0) = (c.x as u32, c.y as u32);
|
|
let mut sum = 0.0;
|
|
for v in v0 - 3..=v0 + 3 {
|
|
for u in u0 - 3..=u0 + 3 {
|
|
sum += (first.dn_at(u, v, 1) as f64 - d.bias) * d.gain;
|
|
}
|
|
}
|
|
sum
|
|
})
|
|
.collect();
|
|
let (m, s) = (mean(&sums), variance(&sums).sqrt());
|
|
let want = signal + 49.0 * dark;
|
|
assert!(
|
|
(m / want - 1.0).abs() < 0.02,
|
|
"{name}: mean {m} e-, expected {want}"
|
|
);
|
|
assert!(
|
|
(s / spread - 1.0).abs() < 0.12,
|
|
"{name}: spread {s} e-, expected {spread} (SNR {})",
|
|
signal / spread
|
|
);
|
|
}
|
|
}
|