From 85cc4273be9d62ba486a2d643b668660f06ee2bf Mon Sep 17 00:00:00 2001 From: dev-null-undefined Date: Sat, 10 Oct 2026 22:29:41 +0200 Subject: [PATCH] LR-EDIT-OPTICS-PROFILE: Sony ARW distortion correction from the embedded table New vendor::arw_lens_rm6 reads the raw IFD's lens tables (0x7037 distortion, 0x7032 vignetting, 0x7035 lateral CA, their flags) and turns the distortion table into a one-plane WarpRectilinear in OpcodeList3, with the conventions measured on 54 ILCE-7RM6 photos: knots at i/15 of the full-frame crop half diagonal, S(r) = r (1 + v(r)/15000), frame-filling zoom, centred on the default crop. Only the ILCE-7RM6 is enabled; vignetting is read, not applied. The radial least-squares fit is split out as vendor::fit_radial_warp (RW2 results unchanged); its error gate can be limited to the radii the crop shows, so a padded active area (lossless files) doesn't refuse a fit over pixels that are never rendered. It sits next to main's vendor::arw_lens (#498), which holds the ILCE-7RM4A conventions measured against Sony's exports (unit 16384, table indexed at the output radius); the two are selected by camera model (disjoint allowlists). arw::decode chains this module's opcodes after arw_lens::distortion, with the same active area and the same crop (the DNG-style default crop when the file has one, which on the ILCE-7RM6 equals Sony's crop tags), and like it only for CFA data (not the downsized linear-RGB lossless sizes, where the knot spacing in full-frame pixels would not hold). The corpus test checks the conversion on the two lossless RawDB files (CC0). The APS-C one is pinned in xtask RAW_SAMPLES with its sha256; the full-frame one is not (at most two pins per PR, and #625 pins two Compressed ARW 6.0 files), so its case runs only where the file is present. RENDER_CACHE_VERSION 62 -> 63. Co-Authored-By: Claude --- crates/engine/src/media.rs | 2 +- crates/raw/src/opcodes.rs | 2 +- crates/raw/src/vendor/arw.rs | 43 ++- crates/raw/src/vendor/arw_lens_rm6.rs | 371 ++++++++++++++++++++++++++ crates/raw/src/vendor/mod.rs | 40 ++- crates/raw/tests/corpus.rs | 63 +++++ xtask/src/main.rs | 7 + 7 files changed, 516 insertions(+), 12 deletions(-) create mode 100644 crates/raw/src/vendor/arw_lens_rm6.rs diff --git a/crates/engine/src/media.rs b/crates/engine/src/media.rs index cf17c3264..5755a2a43 100644 --- a/crates/engine/src/media.rs +++ b/crates/engine/src/media.rs @@ -62,7 +62,7 @@ impl SettingsHashes { } /// Bump when the pipeline's output changes, to invalidate cached thumbnails. -pub const RENDER_CACHE_VERSION: u64 = 62; +pub const RENDER_CACHE_VERSION: u64 = 63; /// Thumbnails render at one of these long edges (so window/cell size changes reuse the cache). pub const THUMB_SIZES: [usize; 4] = [128, 256, 384, 512]; diff --git a/crates/raw/src/opcodes.rs b/crates/raw/src/opcodes.rs index 1d3a2c546..483ac5c4a 100644 --- a/crates/raw/src/opcodes.rs +++ b/crates/raw/src/opcodes.rs @@ -454,7 +454,7 @@ fn gain_at(points_v: u32, points_h: u32, spacing: [f64; 2], origin: [f64; 2], ma } /// Radius normaliser: distance from the centre to the farthest image corner. -fn max_radius(cx: f64, cy: f64, w: f64, h: f64) -> f64 { +pub(crate) fn max_radius(cx: f64, cy: f64, w: f64, h: f64) -> f64 { [(0.0, 0.0), (w, 0.0), (0.0, h), (w, h)].iter().map(|&(x, y)| (x - cx).hypot(y - cy)).fold(0.0, f64::max).max(1e-9) } diff --git a/crates/raw/src/vendor/arw.rs b/crates/raw/src/vendor/arw.rs index 59a9cf59e..5a23218c0 100644 --- a/crates/raw/src/vendor/arw.rs +++ b/crates/raw/src/vendor/arw.rs @@ -2,7 +2,8 @@ //! //! Sources: H. Dietz, "Sony ARW2 Compression: Artifacts And Credible Repair" (IS&T Electronic Imaging 2016) for //! the ARW2 ("cRAW") scheme, and the ExifTool Sony tag-name documentation for the meaning of the raw-IFD tags -//! (`0x7010` tone-curve thresholds, `0x7310` black levels, `0x7313` WB levels, `0x74c7/0x74c8` crop). +//! (`0x7010` tone-curve thresholds, `0x7310` black levels, `0x7313` WB levels, `0x74c7/0x74c8` crop, `0x7031`–`0x7037` +//! lens-correction tables: see [`super::arw_lens`] for the ILCE-7RM4A and [`super::arw_lens_rm6`] for the ILCE-7RM6). //! //! ARW2, as described in the paper: //! 1. Sensor values are tone-mapped to 11-bit codes by a five-segment piecewise-linear curve whose step doubles at @@ -770,7 +771,9 @@ pub(crate) fn decode(bytes: &[u8], mode: Mode) -> Result { // The distortion table is centred on the DNG-style default crop when the file has one: on the ILCE-7RM4A // it starts at x = 32 where Sony's crop tags (the image crop, unchanged) start at 0, and the warp was // validated against Sony's exports with that centre. The opcode's centre is in active-area coordinates, so - // it stays on the optical centre whichever crop frames the image. + // it stays on the optical centre whichever crop frames the image. The ILCE-7RM6's tables follow other + // conventions, measured against its camera JPEGs ([`super::arw_lens_rm6`]); the two allowlists are + // disjoint, so at most one of them produces an opcode. opcodes: OpcodeLists { // DefaultCrop can retain native geometry while Sony's tags select an aspect crop. // Both the image framing and the table geometry must stay within the validated 3:2 scope. @@ -778,7 +781,11 @@ pub(crate) fn decode(bytes: &[u8], mode: Mode) -> Result { Vec::new() } else { let geometry = dng_default_crop(raw, w, h).unwrap_or(crop); - super::arw_lens::distortion(&model, raw, Rect::new(0, 0, w, h), geometry).into_iter().collect() + let active = Rect::new(0, 0, w, h); + super::arw_lens::distortion(&model, raw, active, geometry) + .into_iter() + .chain(super::arw_lens_rm6::lens_opcodes(raw, &model, active, geometry)) + .collect() }, ..Default::default() }, @@ -867,6 +874,36 @@ mod tests { TiffWriter::default().write(&[raw]).unwrap() } + /// The ILCE-7RM6's distortion table reaches `OpcodeList3`, in the headers-only probe and in the full decode; + /// another body's table doesn't (see [`super::super::arw_lens_rm6`]). + #[test] + fn lens_tables_reach_opcode_list3() { + use lightcraft_tiff::{IfdBuilder, ImageData, TiffWriter, Value}; + let file = |model: &str| { + let mut raw = IfdBuilder::new(); + raw.set(t::MAKE, Value::Ascii("SONY".into())); + raw.set(t::MODEL, Value::Ascii(model.into())); + // 3:2, the aspect the opcode is limited to + raw.set(t::IMAGE_WIDTH, Value::Long(vec![96])); + raw.set(t::IMAGE_LENGTH, Value::Long(vec![64])); + raw.set(t::BITS_PER_SAMPLE, Value::Short(vec![12])); + raw.set(t::SAMPLES_PER_PIXEL, Value::Short(vec![1])); + raw.set(t::PHOTOMETRIC, Value::Short(vec![photometric::CFA])); + raw.set(t::COMPRESSION, Value::Short(vec![32767])); + raw.set(TONE_CURVE, Value::Short(vec![8000, 10400, 12900, 14100])); + let table = vec![16, 0, 2, 7, 16, 28, 43, 61, 83, 109, 139, 175, 217, 267, 325, 395, 479]; + raw.set(super::super::arw_lens_rm6::DISTORTION_PARAMS, Value::SShort(table)); + raw.set_image(ImageData::Strips { rows_per_strip: 64, strips: vec![vec![0u8; 96 * 64]] }); + TiffWriter::default().write(&[raw]).unwrap() + }; + for mode in [Mode::Header, Mode::Full] { + let list3 = |model| decode(&file(model), mode).unwrap().opcodes.list3; + let ops = list3("ILCE-7RM6"); + assert!(matches!(ops.as_slice(), [crate::Opcode::WarpRectilinear { planes, .. }] if planes.len() == 1), "{mode:?}: {ops:?}"); + assert!(list3("ILCE-7M4").is_empty(), "{mode:?}"); + } + } + #[test] fn strip_larger_or_smaller_than_one_byte_per_pixel_is_the_packed_variant() { for strip in [32 * 4 + 8, 32 * 4 * 5 / 4 - 1, 32 * 4 - 1] { diff --git a/crates/raw/src/vendor/arw_lens_rm6.rs b/crates/raw/src/vendor/arw_lens_rm6.rs new file mode 100644 index 000000000..ba9aaaa9d --- /dev/null +++ b/crates/raw/src/vendor/arw_lens_rm6.rs @@ -0,0 +1,371 @@ +//! Sony ARW lens corrections from the tables in the raw IFD, as the equivalent DNG `OpcodeList3` opcodes, so they +//! reach the same "profile corrections" (`optics.lens_profile`) as a DNG's embedded warp. +//! +//! These are the ILCE-7RM6 conventions, measured against the camera JPEG. [`super::arw_lens`] holds the ILCE-7RM4A +//! conventions, measured against Sony's Imaging Edge exports (unit 16384, the table indexed at the output radius). +//! Each is selected by camera model; no single unit fits every lens under either. See both docs: +//! docs/sony-ilce7rm6-lens-corrections.md (→ Two conventions) and docs/sony-lens-corrections.md. +//! +//! Sony records three count-prefixed signed 16-bit tables in the raw IFD, each with an in-camera on/off tag (ExifTool +//! Sony tag names): distortion [`DISTORTION_PARAMS`] (flag [`DISTORTION_FLAG`]), vignetting [`VIGNETTING_PARAMS`] +//! ([`VIGNETTING_FLAG`]) and lateral chromatic aberration [`CA_PARAMS`] ([`CA_FLAG`]). The first value is the knot +//! count: 16 / 16 / 32 in full-frame files of the ILCE-7RM6, 11 / 11 / 22 in its APS-C crop mode, the rest padding; +//! the CA table holds red's knots, then blue's. +//! +//! Their meaning was measured by registering blocks of the camera JPEG in each file on LightCraft's uncorrected render +//! (and raw data, for CA) over 54 ILCE-7RM6 photos with four lenses, 41 of them shot with in-camera distortion +//! correction on; nothing else (clean room). The conventions, with their worst residuals at a 1440 px long edge: +//! - knot `i` sits at sensor radius `i · s` from the default-crop centre, `s` = the body's full-frame default-crop +//! half diagonal / 15 (the same pixel spacing in APS-C mode: the knots are fixed on the sensor); linear in between, +//! the last segment extended beyond the last knot. +//! - distortion: a corrected radius `r` (px) reads the source radius `S(r) = r · (1 + v(r) / 15000)` (positive: +//! pincushion), and the camera then zooms (output radius = `z` · corrected radius) so the output fills the default +//! crop and reads only inside it: the largest `S(ρ/z)/ρ` over the crop border radii `ρ` is 1. 0.596 px over the 34 +//! measurable corrected photos (budget 0.75 px; corners and edge midpoints 0.570 px). +//! - `WarpRectilinear`'s 4-term odd polynomial fits that map within 2.35·10⁻⁴ of the DNG radius (0.220 px) on every +//! photo, at its 65 fit samples. A table it fits worse than 0.25 px at the samples the crop shows is refused (the +//! photo stays uncorrected rather than wrongly warped). Between samples the error is somewhat larger: 0.071 against +//! 0.049 px for the full-frame 16–35 mm example in the tests. +//! - vignetting: read but never applied. The camera JPEGs keep the uncorrected render's radial falloff (median +//! corner/centre ratio 0.996 over 39 photos) and no radial gain model fits the tables. +//! - lateral CA: red and blue relative to green on the same knots (unit about 1/3.36·10⁶; 0.399 px, borderline); +//! not applied yet. +//! +//! The correction is produced whatever the in-camera flag says. Only bodies listed in [`FULL_FRAME_CROP`] get it: the +//! unit and the knot spacing are verified on the ILCE-7RM6 alone. Another body is added by repeating the measurement +//! on its photos shot with distortion correction on and listing its full-frame default crop. The method and the +//! numbers: docs/sony-ilce7rm6-lens-corrections.md. + +use super::fit_radial_warp; +use crate::opcodes::max_radius; +use crate::{Opcode, Rect}; +use lightcraft_tiff::Ifd; +use std::ops::RangeInclusive; + +pub(crate) const DISTORTION_PARAMS: u16 = 0x7037; +pub(crate) const VIGNETTING_PARAMS: u16 = 0x7032; +pub(crate) const CA_PARAMS: u16 = 0x7035; +pub(crate) const DISTORTION_FLAG: u16 = 0x7036; +pub(crate) const VIGNETTING_FLAG: u16 = 0x7031; +pub(crate) const CA_FLAG: u16 = 0x7034; + +/// The bodies whose tables are verified (Exif `Model`), with their full-frame default crop (width, height) in px. +const FULL_FRAME_CROP: &[(&str, (f64, f64))] = &[("ILCE-7RM6", (9984.0, 6656.0))]; +/// The distortion table's value of 1 (relative radius); a global least-squares fit over the 34 photos gives 15000.7. +const DISTORTION_UNIT: f64 = 15000.0; +/// Crop border radii scanned by the frame-filling zoom. +const ZOOM_SCAN: usize = 4096; +/// Largest polynomial fit error accepted at the fit samples the crop shows, px at a 1440 px long edge. +const FIT_ERROR_PX_AT_1440: f64 = 0.25; + +/// A file's lens tables: the first `count` values of each (CA: red, blue); `None` when absent or malformed. +/// `vignetting` and `flags` are parsed for the tests and the docs only: don't build an opcode from them (the camera +/// applies no vignetting correction to its JPEGs, and the flags don't decide anything). +#[derive(Debug)] +// only the distortion table becomes an opcode so far +#[cfg_attr(not(test), allow(dead_code))] +pub(crate) struct LensTables { + pub distortion: Option>, + pub vignetting: Option>, + pub ca: Option<(Vec, Vec)>, + /// In-camera on/off of distortion, vignetting and CA correction, as recorded. + pub flags: [Option; 3], +} + +/// The values of a count-prefixed integer table whose count is in `counts` and fits the array. +fn table(raw: &Ifd, tag: u16, counts: RangeInclusive) -> Option> { + let v = raw.value(tag)?; + let count = usize::try_from(v.get_i64(0)?).ok().filter(|n| counts.contains(n) && *n < v.count())?; + (1..=count).map(|i| v.get_i64(i).map(|x| x as f64)).collect() +} + +/// The lens tables of the raw IFD `raw`. +pub(crate) fn read_tables(raw: &Ifd) -> LensTables { + let ca = table(raw, CA_PARAMS, 4..=32).filter(|v| v.len().is_multiple_of(2)).map(|v| { + let (red, blue) = v.split_at(v.len() / 2); + (red.to_vec(), blue.to_vec()) + }); + LensTables { + distortion: table(raw, DISTORTION_PARAMS, 2..=16), + vignetting: table(raw, VIGNETTING_PARAMS, 2..=16), + ca, + flags: [raw.i64(DISTORTION_FLAG), raw.i64(VIGNETTING_FLAG), raw.i64(CA_FLAG)], + } +} + +/// The knot spacing (px) of a verified body: its full-frame default-crop half diagonal / 15. +fn knot_spacing(model: &str) -> Option { + let model = model.trim(); + FULL_FRAME_CROP.iter().find(|(m, _)| *m == model).map(|&(_, (w, h))| w.hypot(h) / 30.0) +} + +/// The knot table at radius `r` px: knot `i` at `i · spacing`, linear in between, the last segment extended beyond +/// the last knot. NaN with fewer than two knots. +fn interpolate(knots: &[f64], spacing: f64, r: f64) -> f64 { + let t = r / spacing; + // `as` saturates, so a huge `t` lands on the last segment; NaN fails `t > 0` + let i = if t > 0.0 { (t as usize).min(knots.len().saturating_sub(2)) } else { 0 }; + match (knots.get(i), knots.get(i + 1)) { + (Some(a), Some(b)) => a + (b - a) * (t - i as f64), + _ => f64::NAN, + } +} + +/// The camera's frame-filling zoom `z` for the source radius map `source` and a crop of half size `x × y`: the +/// largest `source(ρ/z)/ρ` over the border radii `ρ` in [min(x, y), hypot(x, y)] is 1, so the output crop reads +/// only inside the source crop and touches it (at the corner for pincushion, the short-edge midpoint for barrel). +/// Bisection over [0.5, 2]; `None` when that range holds no solution or a ratio is not finite. +fn fill_zoom(source: &dyn Fn(f64) -> f64, x: f64, y: f64) -> Option { + let (near, far) = (x.min(y), x.hypot(y)); + let overshoots = |z: f64| -> Option { + let mut worst = f64::NEG_INFINITY; + for k in 0..=ZOOM_SCAN { + let rb = near + k as f64 * (far - near) / ZOOM_SCAN as f64; + let ratio = source(rb / z) / rb; + if !ratio.is_finite() { + return None; + } + worst = worst.max(ratio); + } + Some(worst > 1.0) + }; + let (mut lo, mut hi) = (0.5, 2.0); + if !overshoots(lo)? || overshoots(hi)? { + return None; + } + for _ in 0..100 { + let z = 0.5 * (lo + hi); + if overshoots(z)? { lo = z } else { hi = z } + } + Some(0.5 * (lo + hi)) +} + +/// The distortion table as a one-plane `WarpRectilinear` relative to the active area `active` (centred on the +/// default crop `crop`, as the camera corrects), with knots `knot_spacing` px apart. `None` when the table is absent +/// or all zero, the zoom or the fit fails, or the fit is worse than 0.25 px at a 1440 px long edge at one of the fit +/// samples within the crop (a padded active area, as in lossless files, reaches radii that are never shown). +pub(crate) fn distortion_opcode(t: &LensTables, knot_spacing: f64, active: Rect, crop: Rect) -> Option { + let v = t.distortion.as_deref().filter(|v| v.iter().any(|&x| x != 0.0))?; + let (aw, ah) = (active.width as f64, active.height as f64); + if aw < 2.0 || ah < 2.0 || crop.width == 0 || crop.height == 0 { + return None; + } + let (cw, ch) = (crop.width as f64, crop.height as f64); + let source = |r: f64| r * (1.0 + interpolate(v, knot_spacing, r) / DISTORTION_UNIT); + let z = fill_zoom(&source, cw / 2.0, ch / 2.0)?; + // DNG's normalisation: pixel-index units, radius 1 = the farthest active-area corner + let (cx, cy) = (crop.x as f64 + (cw - 1.0) / 2.0, crop.y as f64 + (ch - 1.0) / 2.0); + let m = max_radius(cx, cy, aw - 1.0, ah - 1.0); + let max_fit_error = FIT_ERROR_PX_AT_1440 * (cw.max(ch) / 1440.0) / m; + let shown = (cw / 2.0).hypot(ch / 2.0) / m; + let [k0, k1, k2, k3] = fit_radial_warp(&|u| source(u * m / z) / m, max_fit_error, shown)?; + Some(Opcode::WarpRectilinear { planes: vec![[k0, k1, k2, k3, 0.0, 0.0]], center: [cx / (aw - 1.0), cy / (ah - 1.0)] }) +} + +/// The `OpcodeList3` lens corrections of a Sony raw IFD for the body `model` (Exif `Model`); empty for a body whose +/// tables are not verified. +pub(crate) fn lens_opcodes(raw: &Ifd, model: &str, active: Rect, crop: Rect) -> Vec { + let Some(spacing) = knot_spacing(model) else { return Vec::new() }; + distortion_opcode(&read_tables(raw), spacing, active, crop).into_iter().collect() +} + +#[cfg(test)] +mod tests { + use super::*; + use lightcraft_tiff::{Entry, Value}; + + /// The FF 16–35 mm example (corpus `arw-sony-ilce7rm6-ff-land-lossless.arw`, docs/sony-ilce7rm6-lens-corrections.md). + const FF_DISTORTION: [i16; 17] = [16, 0, 2, 7, 16, 28, 43, 61, 83, 109, 139, 175, 217, 267, 325, 395, 479]; + const FF_VIGNETTING: [i16; 17] = [16, 0, 80, 336, 672, 1120, 1648, 2240, 2880, 3600, 4336, 5104, 5872, 6640, 7392, 8160, 8976]; + const FF_CA: [i16; 33] = [ + 32, 242, 234, 224, 212, 196, 176, 154, 132, 108, 86, 66, 46, 30, 16, 8, 4, -276, -258, -238, -212, -184, -148, -114, -74, -40, -6, 20, 42, + 56, 60, 54, 32, + ]; + /// The APS-C example (corpus `arw6-sony-ilce7rm6-apsc-land-hq.arw`): counts 11 / 11 / 22, then padding. + const APSC_DISTORTION: [i16; 17] = [11, 0, -10, -36, -78, -134, -202, -277, -356, -435, -512, -583, -583, -583, -583, -583, -583]; + const APSC_VIGNETTING: [i16; 17] = [11, 0, 48, 208, 496, 880, 1360, 1952, 2624, 3568, 4832, 6176, 0, 0, 0, 0, 0]; + const APSC_CA: [i16; 33] = [ + 22, 790, 764, 734, 702, 666, 624, 582, 544, 510, 484, 470, -432, -404, -372, -330, -286, -230, -174, -114, -58, -6, 34, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, + ]; + /// ARW6 geometry (the decoder's raw image size and default crop). + const FF_ACTIVE: Rect = Rect { x: 0, y: 0, width: 10016, height: 6672 }; + const FF_CROP: Rect = Rect { x: 12, y: 8, width: 9984, height: 6656 }; + const APSC_ACTIVE: Rect = Rect { x: 0, y: 0, width: 6592, height: 4372 }; + const APSC_CROP: Rect = Rect { x: 16, y: 10, width: 6528, height: 4352 }; + + fn knots() -> f64 { + knot_spacing("ILCE-7RM6").unwrap() + } + + fn ifd(tags: Vec<(u16, Value)>) -> Ifd { + let mut entries: Vec = tags.into_iter().map(|(tag, value)| Entry { tag, value, offset: 0 }).collect(); + entries.sort_by_key(|e| e.tag); + Ifd { entries, ..Default::default() } + } + + fn sony_lens_ifd(distortion: &[i16], vignetting: &[i16], ca: &[i16], flags: [u16; 3]) -> Ifd { + ifd(vec![ + (DISTORTION_PARAMS, Value::SShort(distortion.to_vec())), + (VIGNETTING_PARAMS, Value::SShort(vignetting.to_vec())), + (CA_PARAMS, Value::SShort(ca.to_vec())), + (DISTORTION_FLAG, Value::Short(vec![flags[0]])), + (VIGNETTING_FLAG, Value::Short(vec![flags[1]])), + (CA_FLAG, Value::Short(vec![flags[2]])), + ]) + } + + fn tables_with_distortion(v: &[f64]) -> LensTables { + LensTables { distortion: Some(v.to_vec()), vignetting: None, ca: None, flags: [None; 3] } + } + + fn values(table: &[i16]) -> Vec { + table[1..=table[0] as usize].iter().map(|&v| v as f64).collect() + } + + /// The warp's one plane and centre against the recipe's numbers (docs/sony-ilce7rm6-lens-corrections.md). + fn assert_warp(op: Option, center: [f64; 2], k: [f64; 4]) { + let Some(Opcode::WarpRectilinear { planes, center: c }) = op else { panic!("no warp: {op:?}") }; + assert_eq!(planes.len(), 1); + assert!((c[0] - center[0]).abs() < 1e-9 && (c[1] - center[1]).abs() < 1e-9, "centre {c:?}, want {center:?}"); + let p = planes[0]; + for (i, want) in k.iter().enumerate() { + assert!((p[i] - want).abs() < 1e-6, "k{i} {}, want {want} ({p:?})", p[i]); + } + assert_eq!((p[4], p[5]), (0.0, 0.0), "radial only"); + } + + #[test] + fn reads_the_three_tables_and_flags() { + let t = read_tables(&sony_lens_ifd(&FF_DISTORTION, &FF_VIGNETTING, &FF_CA, [0, 2, 1])); + let d = t.distortion.unwrap(); + assert_eq!((d.len(), d[0], d[15]), (16, 0.0, 479.0)); + let v = t.vignetting.unwrap(); + assert_eq!((v.len(), v[1], v[15]), (16, 80.0, 8976.0)); + let (red, blue) = t.ca.unwrap(); + assert_eq!((red.len(), red[0], red[15]), (16, 242.0, 4.0)); + assert_eq!((blue.len(), blue[0], blue[15]), (16, -276.0, 32.0)); + assert_eq!(t.flags, [Some(0), Some(2), Some(1)]); + } + + #[test] + fn apsc_tables_use_only_the_first_count_values() { + let t = read_tables(&sony_lens_ifd(&APSC_DISTORTION, &APSC_VIGNETTING, &APSC_CA, [0, 2, 1])); + let d = t.distortion.unwrap(); + assert_eq!((d.len(), d[10]), (11, -583.0)); + assert_eq!(t.vignetting.unwrap().len(), 11); + let (red, blue) = t.ca.unwrap(); + assert_eq!((red.len(), red[10]), (11, 470.0)); + assert_eq!((blue.len(), blue[0], blue[10]), (11, -432.0, 34.0)); + } + + #[test] + fn missing_or_malformed_tables_yield_no_opcode() { + let mut short = FF_DISTORTION.to_vec(); + short.pop(); + let mut count_17 = FF_DISTORTION.to_vec(); + count_17[0] = 17; + count_17.push(500); + let mut too_large = FF_DISTORTION.to_vec(); + too_large[0] = 600; + for (what, table) in [ + ("count 15 with 3 values", vec![15, 0, 2, 7]), + ("count 17", count_17), + ("empty", vec![]), + ("count 16 with 15 values", short), + ("count 1", vec![1, 5]), + ("count 0", vec![0, 5, 6]), + ("negative count", vec![-16, 0, 2, 7]), + ("count larger than the array", too_large), + ] { + let ifd = ifd(vec![(DISTORTION_PARAMS, Value::SShort(table))]); + let t = read_tables(&ifd); + assert!(t.distortion.is_none(), "{what}"); + assert!(distortion_opcode(&t, knots(), FF_ACTIVE, FF_CROP).is_none(), "{what}"); + assert!(lens_opcodes(&ifd, "ILCE-7RM6", FF_ACTIVE, FF_CROP).is_empty(), "{what}"); + } + // not an integer array, or absent + for value in [Value::Ascii("16 0 2 7".into()), Value::Double(vec![2.0, 0.0, 1.0])] { + assert!(read_tables(&ifd(vec![(DISTORTION_PARAMS, value)])).distortion.is_none()); + } + let none = read_tables(&ifd(vec![])); + assert!(none.distortion.is_none() && none.vignetting.is_none() && none.ca.is_none() && none.flags == [None; 3]); + // the CA table needs an even count of 4..=32 + for ca in [vec![3, 1, 2, 3], vec![34; 35], vec![2, 1, 2]] { + assert!(read_tables(&ifd(vec![(CA_PARAMS, Value::SShort(ca))])).ca.is_none()); + } + } + + #[test] + fn hostile_tables_are_rejected() { + let folded = [0.0, 2000.0, 1000.0, 3000.0, 2000.0, 4000.0, 3000.0, 5000.0, 4000.0, 6000.0, 5000.0, 7000.0, 6000.0, 8000.0, 7000.0, 9000.0]; + for (what, table) in [ + ("huge", vec![(i64::MAX / 4) as f64; 16]), + ("most negative", vec![i64::MIN as f64; 16]), + ("zig-zag", folded.to_vec()), + ("all zero", vec![0.0; 16]), + ("one knot", vec![100.0]), + ("no knots", vec![]), + ("not finite", vec![0.0, f64::NAN, 20.0]), + ("infinite", vec![0.0, f64::INFINITY, 20.0]), + ] { + assert!(distortion_opcode(&tables_with_distortion(&table), knots(), FF_ACTIVE, FF_CROP).is_none(), "{what}"); + } + // the same through the IFD, as 64-bit values + for v in [i64::MAX / 4, i64::MIN, i64::MAX] { + let mut table = vec![v; 17]; + table[0] = 16; + let ifd = ifd(vec![(DISTORTION_PARAMS, Value::SLong8(table))]); + assert!(lens_opcodes(&ifd, "ILCE-7RM6", FF_ACTIVE, FF_CROP).is_empty(), "{v}"); + } + // degenerate geometry + let ff = tables_with_distortion(&values(&FF_DISTORTION)); + for (what, knots, active, crop) in [ + ("empty crop", knots(), FF_ACTIVE, Rect::new(12, 8, 0, 6656)), + ("1 × 1 active area", knots(), Rect::new(0, 0, 1, 1), Rect::new(0, 0, 1, 1)), + ("knot spacing NaN", f64::NAN, FF_ACTIVE, FF_CROP), + ("knot spacing 0", 0.0, FF_ACTIVE, FF_CROP), + ] { + assert!(distortion_opcode(&ff, knots, active, crop).is_none(), "{what}"); + } + } + + #[test] + fn ff_example_table_becomes_the_measured_warp() { + let t = tables_with_distortion(&values(&FF_DISTORTION)); + assert_warp( + distortion_opcode(&t, knots(), FF_ACTIVE, FF_CROP), + [0.4996005991, 0.5], + [0.9714495823, 0.022508032, 0.0006090085237, 0.005683284548], + ); + } + + #[test] + fn apsc_example_table_becomes_the_measured_warp() { + let t = tables_with_distortion(&values(&APSC_DISTORTION)); + let op = distortion_opcode(&t, knots(), APSC_ACTIVE, APSC_CROP); + assert_warp(op, [0.4975724473, 0.5], [1.01627551, -0.0618606692, 0.02768349621, -0.005595798468]); + } + + /// The lossless files' padded active area (10240 × 7168) puts the farthest active-area corner (DNG radius 1) well + /// beyond the crop's (u 0.925): the gate looks only at the radii the crop shows. The FF example × 1.5 fits within + /// 0.12 px at 1440 at the samples there, 0.27 px at all 65. + #[test] + fn padded_active_area_gates_only_the_displayed_radii() { + let strong: Vec = values(&FF_DISTORTION).iter().map(|v| v * 1.5).collect(); + let padded = Rect { x: 0, y: 0, width: 10240, height: 7168 }; + let op = distortion_opcode(&tables_with_distortion(&strong), knots(), padded, FF_CROP); + assert!(matches!(op, Some(Opcode::WarpRectilinear { .. })), "{op:?}"); + } + + #[test] + fn unknown_model_gets_no_opcodes() { + let ifd = sony_lens_ifd(&FF_DISTORTION, &FF_VIGNETTING, &FF_CA, [0, 2, 1]); + assert!(lens_opcodes(&ifd, "ILCE-7M4", FF_ACTIVE, FF_CROP).is_empty()); + assert!(lens_opcodes(&ifd, "", FF_ACTIVE, FF_CROP).is_empty()); + // the camera's in-camera flag (0 = off here) doesn't matter: the table is applied + let ops = lens_opcodes(&ifd, " ILCE-7RM6 ", FF_ACTIVE, FF_CROP); + assert_eq!(ops, distortion_opcode(&read_tables(&ifd), knots(), FF_ACTIVE, FF_CROP).into_iter().collect::>()); + assert_eq!(ops.len(), 1); + } +} diff --git a/crates/raw/src/vendor/mod.rs b/crates/raw/src/vendor/mod.rs index e5218f390..5512a1169 100644 --- a/crates/raw/src/vendor/mod.rs +++ b/crates/raw/src/vendor/mod.rs @@ -6,6 +6,7 @@ mod arw6_curve; #[cfg(test)] pub(crate) mod arw6_testenc; mod arw_lens; +pub(crate) mod arw_lens_rm6; pub mod cr2; pub mod cr3; pub mod crx; @@ -111,7 +112,7 @@ pub(crate) fn white_from_data(data: &[u16], bits: u32) -> f32 { } } -/// Output radii sampled to fit the inverse mapping. +/// Output radii sampled to fit a radial warp. const DISTORTION_SAMPLES: usize = 64; /// Least-squares radial warp from a camera's forward distortion map. @@ -121,10 +122,10 @@ const DISTORTION_SAMPLES: usize = 64; /// from the centre to the **farthest** corner = 1) so the result is directly a `WarpRectilinear` plane. The map is /// scanned at 4096 points up to a source radius of `3 * t_max`: it must reach `t_max` within that range and rise at /// every scan point (monotonicity is checked only there), otherwise `None`; a NaN at a scan point is rejected the -/// same way, other non-finite values are caught by the fit-error gate. The map is inverted by bisection and -/// `source = u·(k0 + k1u² + k2u⁴ + k3u⁶)` is fitted over output radii `u·t_max`, `u` in [0, 1]: the returned -/// `[k0, k1, k2, k3]` give the source radius in the caller's unit as a function of `u`. `None` too when the worst -/// fit error exceeds `max_fit_error`. +/// same way, other non-finite values by the fit. The map is inverted by bisection and +/// `source = u·(k0 + k1u² + k2u⁴ + k3u⁶)` is fitted over output radii `u·t_max`, `u` in [0, 1] ([`fit_radial_warp`]): +/// the returned `[k0, k1, k2, k3]` give the source radius in the caller's unit as a function of `u`. `None` too when +/// the worst fit error exceeds `max_fit_error`. pub(crate) fn warp_from_forward_map(corrected: &dyn Fn(f64) -> f64, t_max: f64, max_fit_error: f64) -> Option<[f64; 4]> { if !(t_max.is_finite() && t_max > 0.0) { return None; @@ -152,12 +153,23 @@ pub(crate) fn warp_from_forward_map(corrected: &dyn Fn(f64) -> f64, t_max: f64, } 0.5 * (lo + hi) }; + fit_radial_warp(&|u| source(u * t_max), max_fit_error, 1.0) +} + +/// Least-squares `source = u·(k0 + k1u² + k2u⁴ + k3u⁶)` through `source_of_u` sampled at `u = i/64`, i = 0..=64 +/// (`u` = output radius, `WarpRectilinear` units, so the result is directly a plane's `[k0, k1, k2, k3]`). `None` when +/// a sample is not finite, the normal equations are singular, or the worst fit error at the samples with +/// `u ≤ gate_until` (the radii the image shows; 1 = all) exceeds `max_fit_error`. +pub(crate) fn fit_radial_warp(source_of_u: &dyn Fn(f64) -> f64, max_fit_error: f64, gate_until: f64) -> Option<[f64; 4]> { let samples: Vec<(f64, f64)> = (0..=DISTORTION_SAMPLES) .map(|i| { let u = i as f64 / DISTORTION_SAMPLES as f64; - (u, source(u * t_max)) + (u, source_of_u(u)) }) .collect(); + if samples.iter().any(|&(_, v)| !v.is_finite()) { + return None; + } let basis = |u: f64| { let u2 = u * u; [u, u * u2, u * u2 * u2, u * u2 * u2 * u2] @@ -175,7 +187,7 @@ pub(crate) fn warp_from_forward_map(corrected: &dyn Fn(f64) -> f64, t_max: f64, } let k = solve4(a, rhs)?; let fit = |u: f64| basis(u).iter().zip(&k).map(|(f, c)| f * c).sum::(); - let worst = samples.iter().map(|&(u, v)| (fit(u) - v).abs()).fold(0.0, f64::max); + let worst = samples.iter().filter(|&&(u, _)| u <= gate_until).map(|&(u, v)| (fit(u) - v).abs()).fold(0.0, f64::max); (worst <= max_fit_error).then_some(k) } @@ -281,4 +293,18 @@ mod tests { assert!(warp_from_forward_map(&dip, 1.0, 1e-2).is_none()); assert!(warp_from_forward_map(&|_| f64::NAN, 1.0, 1e-3).is_none()); } + + #[test] + fn fit_radial_warp_identity_and_gate() { + let k = fit_radial_warp(&|u| u, 1e-9, 1.0).unwrap(); + assert!((k[0] - 1.0).abs() < 1e-9 && k[1].abs() < 1e-9 && k[2].abs() < 1e-9 && k[3].abs() < 1e-9, "{k:?}"); + assert!(fit_radial_warp(&|u| if u > 0.5 { f64::NAN } else { u }, 1.0, 1.0).is_none(), "a NaN sample"); + assert!(fit_radial_warp(&|u| if u > 0.5 { f64::INFINITY } else { u }, 1.0, 1.0).is_none(), "an infinite sample"); + assert!(fit_radial_warp(&|u: f64| u.sqrt(), 1.0, 1.0).is_some()); + assert!(fit_radial_warp(&|u: f64| u.sqrt(), 0.0, 1.0).is_none(), "not a polynomial: the gate rejects it"); + // a step at u = 1 alone: worst fit error 0.057 over all samples, 0.006 over u <= 0.9 (the fit is the same) + let step = |u: f64| if u > 0.99 { 1.1 } else { u }; + assert!(fit_radial_warp(&step, 0.01, 1.0).is_none()); + assert_eq!(fit_radial_warp(&step, 0.01, 0.9), fit_radial_warp(&step, 1.0, 1.0), "gated only where shown"); + } } diff --git a/crates/raw/tests/corpus.rs b/crates/raw/tests/corpus.rs index 5ca4777b7..60c81fd1f 100644 --- a/crates/raw/tests/corpus.rs +++ b/crates/raw/tests/corpus.rs @@ -626,6 +626,69 @@ fn corpus_sony_16x9_crops() { eprintln!("Sony 16:9 crops checked on {seen} files"); } +/// Sony ILCE-7RM6: the distortion table in the raw IFD (`0x7037`) becomes a `WarpRectilinear` in `OpcodeList3` +/// (`vendor::arw_lens_rm6`). Expected: the conversion recipe measured on 54 ILCE-7RM6 photos, evaluated in f64 for each +/// file's decoded active area (the padded lossless tiles) and default crop. Both files carry distortion correction +/// off in camera; the table is applied anyway. +#[test] +fn corpus_sony_lens_tables_become_opcodes() { + let dir = corpus_root().join("raw"); + // (file, active area, default crop x, y, width, height, opcode centre, G-plane k0..k3) + let cases = [ + // not pinned in `cargo xtask corpus` (at most two pins per PR; #625 pins two Compressed ARW 6.0 files): checked + // only where the file is present (RawDB `ILCE-7RM6_FF_land_RAW.ARW`, CC0) + ( + "arw-sony-ilce7rm6-ff-land-lossless.arw", + (10240, 7168), + (12, 8, 9984, 6656), + [0.488670768629749, 0.4653969582810102], + [0.9717033132764651, 0.023249297636569695, 0.008638979168841135, 0.003091041831374855], + ), + // pinned in `cargo xtask corpus` (RawDB `ILCE-7RM6_APSC_land_RAW.ARW`, CC0) + ( + "arw-sony-ilce7rm6-apsc-land-lossless.arw", + (6656, 4608), + (16, 10, 6528, 4352), + [0.49278737791134486, 0.4743868026915563], + [1.016366854640611, -0.0688824189782517, 0.03615745253291194, -0.009438825010476538], + ), + ]; + let mut seen = 0; + for (name, (aw, ah), crop, center, k) in cases { + let bytes = match std::fs::read(dir.join(name)) { + Ok(bytes) => bytes, + Err(e) if e.kind() == std::io::ErrorKind::NotFound => { + eprintln!("skip: {name} absent"); + continue; + } + Err(e) => panic!("{name}: {e}"), + }; + let img = decode(&bytes).unwrap_or_else(|e| panic!("{name}: {e}")); + let (a, c) = (img.active_area, img.crop); + eprintln!("{name}: active area {a:?}, crop {c:?}"); + assert_eq!((a.x, a.y, a.width, a.height), (0, 0, aw, ah), "{name}: active area"); + assert_eq!((c.x, c.y, c.width, c.height), crop, "{name}: default crop"); + assert_eq!(probe_info(&bytes).unwrap().opcodes, img.opcodes, "{name}: the headers-only probe reads the same correction"); + let warps: Vec<_> = img + .opcodes + .list3 + .iter() + .filter_map(|op| match op { + lightcraft_raw::Opcode::WarpRectilinear { planes, center } => Some((planes, center)), + _ => None, + }) + .collect(); + let [(planes, got)] = warps.as_slice() else { panic!("{name}: {:?}", img.opcodes.list3) }; + assert_eq!(planes.len(), 1, "{name}: one plane (G)"); + assert!((got[0] - center[0]).abs() < 1e-6 && (got[1] - center[1]).abs() < 1e-6, "{name}: centre {got:?}, want {center:?}"); + for (i, want) in k.iter().enumerate() { + assert!((planes[0][i] - want).abs() < 1e-6, "{name}: k{i} {}, want {want}", planes[0][i]); + } + seen += 1; + } + eprintln!("Sony lens tables checked on {seen} files"); +} + /// The embedded JPEG as linear RGB, reduced to `gw × gh` cells. fn jpeg_cells(jpeg: &[u8], gw: usize, gh: usize) -> Vec<[f64; 3]> { use zune_core::{bytestream::ZCursor, colorspace::ColorSpace, options::DecoderOptions}; diff --git a/xtask/src/main.rs b/xtask/src/main.rs index 08ef435c5..28fb3a38b 100644 --- a/xtask/src/main.rs +++ b/xtask/src/main.rs @@ -448,6 +448,13 @@ const RAW_SAMPLES: &[(&str, &str, &str)] = &[ "https://raw.pixls.us/data/Sony/DSLR-A900/_DSC7969.ARW", "ac7c1532df77c321e8010aa1be60c1b9f245b1ac0db6f38fe617a73aae81af50", ), + // Sony ILCE-7RM6 lossless, CC0 1.0 (RawDB): the APS-C file, for its lens tables (`vendor::arw_lens_rm6`). The + // full-frame lossless file is not pinned (at most two pins per PR; #625 pins two Compressed ARW 6.0 files). + ( + "arw-sony-ilce7rm6-apsc-land-lossless.arw", + "https://rawdb.dnglab.org/api/download/Sony/ILCE-7RM6/raw_modes/ILCE-7RM6_APSC_land_RAW.ARW", + "7aaf0b3941509d7364a6f03844783155ebc912d246d1c456e5cdb4693ec09fbc", + ), // CR2 colour-filter layouts differ by model (issue #85): CR2CFAPattern 3 (GBRG) and 1 (RGGB) samples ( "cr2-canon-40d.cr2",