diff --git a/Cargo.lock b/Cargo.lock index 5be8eae71a..151a5eb0bf 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -2031,7 +2031,10 @@ name = "filmcraft-dng" version = "0.16.0" dependencies = [ "filmcraft-dng", + "image", "rayon", + "serde", + "serde_json", "thiserror 2.0.21", ] diff --git a/crates/dng/Cargo.toml b/crates/dng/Cargo.toml index 1fce4930ce..ca449ff432 100644 --- a/crates/dng/Cargo.toml +++ b/crates/dng/Cargo.toml @@ -11,13 +11,20 @@ default = ["threads"] threads = ["dep:rayon"] # A small DNG and lossless JPEG writer for synthesising test files (not used by the app). writer = [] +# serde derives on the raw settings (`ProcessOptions`) for UI state and project files. +serde = ["dep:serde"] [dependencies] rayon = { workspace = true, optional = true } +# Baseline (lossy) JPEG for Compression 34892 tiles: the JPEG decoder the workspace already uses +# (crates/codecs, crates/media), pure Rust. +image = { version = "0.25", default-features = false, features = ["jpeg"] } +serde = { workspace = true, optional = true } thiserror.workspace = true [dev-dependencies] -filmcraft-dng = { path = ".", features = ["writer"] } +filmcraft-dng = { path = ".", features = ["writer", "serde"] } +serde_json.workspace = true [lints] workspace = true diff --git a/crates/dng/README.md b/crates/dng/README.md index c8678c1d37..964ba345cc 100644 --- a/crates/dng/README.md +++ b/crates/dng/README.md @@ -1,7 +1,8 @@ # filmcraft-dng -Clean-room DNG and CinemaDNG camera RAW decoder in pure Rust. Layer L0: depends only on -`thiserror` and (optional, feature `threads`, default on) `rayon`. Builds for +Clean-room DNG and CinemaDNG camera RAW decoder in pure Rust. Layer L0: depends on `thiserror`, +`image` (its baseline JPEG decoder only, for lossy JPEG DNG; the workspace already uses it in +`codecs` and `media`) and (optional, feature `threads`, default on) `rayon`. Builds for `wasm32-unknown-unknown`; no `unsafe`. Implemented from the public specifications below only; the source of no other raw decoder (dcraw, LibRaw, rawspeed, darktable, RawTherapee, FFmpeg, the Adobe DNG SDK…) was read. LibRaw's `unprocessed_raw` sample program is run as an external test oracle @@ -11,17 +12,18 @@ when it is installed, never linked. | Document | Edition | Used for | |---|---|---| -| Adobe Digital Negative (DNG) Specification | 1.4.0.0 (2012) – 1.7.1.0 (2023) | raw IFD selection (NewSubFileType 0), CFA / LinearRaw data, BitsPerSample packing, LinearizationTable, BlackLevel (+ RepeatDim, DeltaH / DeltaV), WhiteLevel, ActiveArea, DefaultCrop, the colour model of chapter 6 (ColorMatrix, CameraCalibration, AnalogBalance, ForwardMatrix, AsShotNeutral / AsShotWhiteXY, CalibrationIlluminant interpolation, BaselineExposure) | +| Adobe Digital Negative (DNG) Specification | 1.4.0.0 (2012) – 1.7.1.0 (2023) | Compression 34892 (lossy JPEG: 8-bit LinearRaw only), raw IFD selection (NewSubFileType 0), CFA / LinearRaw data, BitsPerSample packing, LinearizationTable, BlackLevel (+ RepeatDim, DeltaH / DeltaV), WhiteLevel, ActiveArea, DefaultCrop, the colour model of chapter 6 (ColorMatrix, CameraCalibration, AnalogBalance, ForwardMatrix, AsShotNeutral / AsShotWhiteXY, CalibrationIlluminant interpolation, BaselineExposure), opcode lists (chapter 7: OpcodeList1–3 and opcodes 1–14) | | TIFF Revision 6.0 | June 1992 | header, IFDs, field types, strips and tiles, PlanarConfiguration, row padding | | TIFF/EP, ISO 12234-2 | 2001 | CFARepeatPatternDim, CFAPattern | | TIFF Technical Note 1 (SubIFDs) | 1993 | the SubIFD tree | | ITU-T T.81 \| ISO/IEC 10918-1 | 09/1992 | lossless JPEG (process 14, SOF3): Annex B syntax, Annex C Huffman tables, F.2.2.3 decoding, Annex H prediction and point transform | -| CinemaDNG Image Data Format Specification | 1.1.0.0 (2009) | frames as numbered DNG files, FrameRate, TimeCodes, ReelName tags | +| CinemaDNG Image Data Format Specification | 1.0.0.0 (September 2009) | frames as numbered DNG files, FrameRate, TimeCodes, ReelName tags | | SMPTE ST 12-1 | 2014 | time code BCD packing (via the CinemaDNG TimeCodes tag) | | CIPA DC-008 (Exif) | 2.3 | LightSource codes → illuminant temperatures | | Malvar, He, Cutler, "High-quality linear interpolation for demosaicing of Bayer-patterned color images", IEEE ICASSP | 2004 | the Bayer demosaic filters | | McCamy, "Correlated color temperature as an explicit function of chromaticity coordinates", Color Res. Appl. 17(2) | 1992 | white xy → temperature for the calibration interpolation | | IEC 61966-2-1 (sRGB) | 1999 | Rec. 709 / sRGB primaries and D65 white | +| Krystek, "An algorithm to calculate correlated colour temperature", Color Res. Appl. 10(1) | 1985 | the Planckian locus in CIE 1960 uv for temperature / tint white balance | ## API @@ -33,9 +35,19 @@ if filmcraft_dng::sniff(&bytes) { let img = dng.develop()?; // linear Rec. 709 RGBA f32 (alpha 1), info.width × info.height } let img = filmcraft_dng::decode(&bytes)?; // parse + develop +// raw settings: white balance (as shot or temperature / tint), exposure, highlights, demosaic +// quality, output space, and a target size for a reduced-resolution fast path +let o = filmcraft_dng::ProcessOptions { exposure: 0.5, ..Default::default() }; +let img = filmcraft_dng::decode_with(&bytes, &o)?; ``` -`writer` (feature `writer`, tests only) writes synthetic DNGs: uncompressed at any depth or +`ProcessOptions::default()` develops exactly like `decode`. `Highlights::Recover` is accepted but +not implemented yet (it develops like `Clip`). `Dng::as_shot_temperature_tint` gives the as-shot +white as kelvin and tint, the starting point of a white-balance control; temperature and tint +are a point on the Planckian locus (Krystek 1985) and a CIE 1960 Δuv (tint 150 = Δuv +0.05, a +greener white, so a more magenta picture; D65 is about tint +10). + +`writer` (feature `writer`, tests only) writes synthetic DNGs: uncompressed at any depth, lossy JPEG (through `image`'s encoder) or lossless JPEG with any predictor, strips or tiles, either byte order, an optional preview IFD with the raw image in a SubIFD, and arbitrary extra tags. @@ -45,7 +57,9 @@ the raw image in a SubIFD, and arbitrary extra tags. cycles end the walk). The main image is the largest NewSubFileType 0 IFD with a CFA or LinearRaw photometric interpretation. - **Raw data:** 1–16-bit unsigned samples, uncompressed (8 and 16-bit in file byte order, other - depths packed most significant bit first) or lossless JPEG (all predictors 1–7, 1–4 components, + depths packed most significant bit first), lossy JPEG (Compression 34892: baseline DCT tiles or + strips of one or three components, 8-bit LinearRaw only, as DNG 1.4 allows; the JPEG's own + colour transform applies) or lossless JPEG (all predictors 1–7, 1–4 components, point transform, restart intervals of whole lines), in strips or tiles; tiles decode in parallel. - **Processing:** ActiveArea, LinearizationTable, per-site black levels with row / column deltas, @@ -53,6 +67,22 @@ the raw image in a SubIFD, and arbitrary extra tags. (Malvar–He–Cutler for 2 × 2 Bayer patterns of any phase, bilinear for other patterns up to 8 × 8), DefaultCrop, the DNG colour model to XYZ D50 and Bradford to linear Rec. 709, and BaselineExposure. LinearRaw (3 samples) skips the demosaic; monochrome LinearRaw is grey. +- **Opcode lists (DNG chapter 7), opt-in (`ProcessOptions::apply_opcodes`, off by default + until the interpretations below are compared with real camera files):** OpcodeList1 on the stored samples (0–65535), OpcodeList2 on + the linear reference values before white balance, OpcodeList3 on the demosaiced camera values + (white balance taken out, 0–1), values clipped after every opcode. WarpRectilinear, + WarpFisheye, WarpRectilinear2 (with the skip rule for a following warp), FixVignetteRadial, + FixBadPixelsConstant / List, TrimBounds, MapTable, MapPolynomial, GainMap, DeltaPerRow / + Column, ScalePerRow / Column. Unknown opcodes, and opcodes from a DNG version above 1.7.1.0, + are skipped when optional and an error naming them when not; preview-only opcodes are skipped + on the reduced-resolution path. Where the specification leaves room: OpcodeList1 uses stored + image coordinates, OpcodeList2 and 3 ActiveArea coordinates; GainMap positions are row / height + and column / width of the whole picture; Delta and MapPolynomial values are in the 0–1 range at + every stage (×65535 on OpcodeList1), MapTable indexes 16-bit values; warps resample with a + Catmull–Rom cubic; warps on a CFA mosaic and TrimBounds before OpcodeList3 are not applied + (an error unless optional). A TrimBounds makes the developed picture smaller than `Info`'s + size. At most 256 opcodes per list; counts are checked against the parameter bytes before + anything is allocated. - **Output:** scene-linear Rec. 709 floats: a white the sensor clipped is 1.0 × 2^BaselineExposure; brighter values are kept (FilmCraft's compositor works in linear float). Negative results of the colour matrix are clipped to 0. Black levels scale per pixel by (white − black). @@ -60,16 +90,36 @@ the raw image in a SubIFD, and arbitrary extra tags. values to XYZ and a Bradford adaptation moves the as-shot white to D50. Decoders in the dcraw tradition (LibRaw included) scale the camera channels instead, so under a white far from D65 their saturated colours differ from ours (greys agree). -- **Speed:** a 4096 × 2160 12-bit frame develops in about 80–100 ms on 16 threads (raw decode - ~15 ms, the rest demosaic and colour). +- **Speed:** a 4096 × 2160 12-bit lossless-JPEG frame decodes and develops about 2–3 times + faster than in the first version (best of 9 on a shared 16-thread machine: tiles 256 × 256 + about 82 ms on 1 thread and 21 ms on 16, one strip 30 ms of serial entropy decoding plus the + processing). Lossless JPEG decodes a code and its difference bits with one 13-bit lookup and + refills its bit buffer eight bytes at a time; tiles decode in parallel straight into place (a + single strip is used without a copy); without black-level deltas a lookup table per CFA / + black-level site replaces the normalisation arithmetic; demosaic, crop and colour run row by + row without intermediate pictures. Output is bit-identical to the first version + (`tests/it/equivalence.rs` compares f32 bit patterns with the kept `reference` module). +- **Reduced resolution:** with `ProcessOptions::target_size`, a Bayer, LinearRaw or monochrome + picture develops at 1/2, 1/4 or 1/8 size (the smallest that still covers the target) by + averaging each colour over 2 × 2 (4 × 4, 8 × 8) blocks of the crop, without a demosaic: about + three times less processing at half size. +- `cargo run --release -p filmcraft-dng --example bench` times the current pipeline against the + reference at 1, 4 and 16 threads, interleaved in one run. ## Not supported (errors name them) -Lossy JPEG and JPEG XL DNG, deflate / floating-point DNG (HDR merges), planar raw data, colour -filters other than red / green / blue (CMY, RGBE), non-rectangular CFA layouts, opcode lists -(warp, vignette, bad pixels: ignored), DefaultScale and BestQualityScale (non-square pixels are -shown square), DNG 1.6 IlluminantData for illuminant 255 (calibration 1 is used), lossless JPEG -restart intervals that split a line, and MXF-wrapped CinemaDNG. +JPEG XL DNG, deflate / floating-point DNG (HDR merges), planar raw data, colour +filters other than red / green / blue (CMY, RGBE), non-rectangular CFA layouts, DefaultScale and +BestQualityScale (non-square pixels are shown square), DNG 1.6 IlluminantData for illuminant 255 +(calibration 1 is used) and lossless JPEG restart intervals that split a line. This crate reads +DNG frames only; MXF-wrapped CinemaDNG is unwrapped by the MXF reader (`crates/mxf`, #1078), +which hands the frames to it. + +Known deviation: the rescale step divides by WhiteLevel minus each pixel's own black level; +DNG chapter 5 uses the plane's largest black level. The two agree unless the black levels +differ across a plane (BlackLevelRepeatDim sites of one colour or BlackLevelDeltaH / V), where +our scale differs by the black-level spread over the white level (well under 1 % for usual +files). Kept for now so the faster pipeline stays bit-identical to the first one. ## Tests @@ -80,6 +130,11 @@ restart intervals that split a line, and MXF-wrapped CinemaDNG. white balance, ForwardMatrix, BaselineExposure, ActiveArea, DefaultCrop, LinearRaw); CinemaDNG FrameRate / TimeCodes / ReelName; hostile files (every truncation, 12 000 mutations under `catch_unwind`, cyclic IFDs, huge sizes, bad offsets, zero denominators). +- `tests/it/lossy.rs`: lossy JPEG tiles and strips with one and three components, samples equal to + the JPEG decoder's, a LinearizationTable expanding 8-bit codes, CFA and 12-bit lossy files + refused, a tiny file claiming a huge picture, truncations and 3000 mutations. +- `tests/it/opcodes.rs`, `tests/it/equivalence.rs`: opcode lists per stage; the faster pipeline + against the reference. - `tests/it/oracle_libraw.rs`: LibRaw's `unprocessed_raw` reads the same samples from our synthetic files (uncompressed 8–16-bit in both byte orders, lossless JPEG predictors 1–7 with one and two components). Skipped when it is not installed. diff --git a/crates/dng/examples/bench.rs b/crates/dng/examples/bench.rs new file mode 100644 index 0000000000..c994a99ec1 --- /dev/null +++ b/crates/dng/examples/bench.rs @@ -0,0 +1,102 @@ +//! Decode + develop speed of a synthetic 4096 × 2160 12-bit lossless-JPEG CinemaDNG frame, the +//! current pipeline against the reference one of PR #1068, interleaved in one run at 1, 4 and 16 +//! threads (rayon pools). `cargo run --release -p filmcraft-dng --features writer --example bench` +//! (`-- --reps N`, `-- --quick` for a 1920 × 1080 frame, `-- --target 1920x1080` for the reduced +//! path). Times are the best of the repetitions, in ms. + +use std::time::Instant; + +use filmcraft_dng::tiff::tags; +use filmcraft_dng::writer::{Compression, Layout, Options, Value, write}; +use filmcraft_dng::{Dng, ProcessOptions}; + +fn frame(w: usize, h: usize, layout: Layout) -> Vec { + // a smooth scene with sensor-like noise, black 256, white 4095 + let mut s = 0x1234_5678u32; + let px: Vec = (0..w * h) + .map(|i| { + s ^= s << 13; + s ^= s >> 17; + s ^= s << 5; + let (x, y) = ((i % w) as f32 / w as f32, (i / w) as f32 / h as f32); + let c = (i % w) % 2 + (i / w) % 2; + let v = 0.15 + 0.6 * (x * 3.0 + y * 2.0 + c as f32).sin().abs() * y; + (256.0 + v * 3500.0 + (s % 48) as f32).min(4095.0) as u16 + }) + .collect(); + let cm = [0.9f64, -0.3, -0.05, -0.4, 1.2, 0.2, -0.05, 0.15, 0.6]; + let o = Options { + width: w, + height: h, + bits: 12, + layout, + compression: Compression::LosslessJpeg { predictor: 1, components: 2, restart_lines: 0, pt: 0 }, + main_tags: vec![ + (tags::COLOR_MATRIX_1, Value::srationals(&cm)), + (tags::CALIBRATION_ILLUMINANT_1, Value::Short(vec![21])), + (tags::AS_SHOT_NEUTRAL, Value::rationals(&[0.55, 1.0, 0.7])), + (tags::FRAME_RATE, Value::SRational(vec![(24, 1)])), + ], + raw_tags: vec![(tags::BLACK_LEVEL, Value::Short(vec![256])), (tags::WHITE_LEVEL, Value::Short(vec![4095]))], + ..Default::default() + }; + write(&px, &o) +} + +/// Best of the repetitions: the least disturbed by other work on a shared machine. +fn best(v: Vec) -> f64 { + v.into_iter().fold(f64::INFINITY, f64::min) +} + +fn ms(t: Instant) -> f64 { + t.elapsed().as_secs_f64() * 1000.0 +} + +fn main() { + let args: Vec = std::env::args().collect(); + let quick = args.iter().any(|a| a == "--quick"); + // profiling: skip the reference + let only_new = args.iter().any(|a| a == "--only-new"); + let reps: usize = args.iter().position(|a| a == "--reps").and_then(|i| args.get(i + 1)).and_then(|v| v.parse().ok()).unwrap_or(7); + let (w, h) = if quick { (1920, 1080) } else { (4096, 2160) }; + let target = args.iter().position(|a| a == "--target").and_then(|i| args.get(i + 1)).and_then(|v| { + let (a, b) = v.split_once('x')?; + Some((a.parse().ok()?, b.parse().ok()?)) + }); + let opts = ProcessOptions { target_size: target, ..Default::default() }; + for (name, layout) in [("tiles 256x256", Layout::Tiles { width: 256, height: 256 }), ("one strip", Layout::Strips { rows: h })] { + let f = frame(w, h, layout); + println!("\n{w}x{h} 12-bit LJ92, {name}, {} MB, options {opts:?}", f.len() / 1_000_000); + println!("loadavg: {}", std::fs::read_to_string("/proc/loadavg").unwrap_or_default().trim()); + println!( + "{:>7} | {:>9} {:>9} {:>9} | {:>9} {:>9} {:>9} | {:>7}", + "threads", "ref raw", "ref dev", "ref all", "new raw", "new dev", "new all", "speedup" + ); + for threads in [1, 4, 16] { + let pool = rayon::ThreadPoolBuilder::new().num_threads(threads).build().unwrap(); + let (mut rr, mut ra, mut nr, mut na) = (vec![], vec![], vec![], vec![]); + pool.install(|| { + for _ in 0..reps { + let d = Dng::parse(&f).unwrap(); + if !only_new { + let t = Instant::now(); + std::hint::black_box(d.reference_raw().unwrap()); + rr.push(ms(t)); + let t = Instant::now(); + std::hint::black_box(d.reference_develop(&opts).unwrap()); + ra.push(ms(t)); + } + let t = Instant::now(); + std::hint::black_box(d.raw().unwrap()); + nr.push(ms(t)); + let t = Instant::now(); + std::hint::black_box(d.develop_with(&opts).unwrap()); + na.push(ms(t)); + } + }); + let (rr, ra, nr, na) = (best(rr), best(ra), best(nr), best(na)); + println!("{threads:>7} | {rr:>9.1} {:>9.1} {ra:>9.1} | {nr:>9.1} {:>9.1} {na:>9.1} | {:>6.2}x", ra - rr, na - nr, ra / na); + } + println!("loadavg: {}", std::fs::read_to_string("/proc/loadavg").unwrap_or_default().trim()); + } +} diff --git a/crates/dng/src/color.rs b/crates/dng/src/color.rs index fb02f2c8d4..d2b9a6e8a0 100644 --- a/crates/dng/src/color.rs +++ b/crates/dng/src/color.rs @@ -101,6 +101,96 @@ pub fn cct((x, y): (f64, f64)) -> f64 { if t.is_finite() { t.clamp(1_000.0, 50_000.0) } else { 50_000.0 } } +/// Tint units per Δuv, with the sign of the usual raw-converter control ("more tint, more +/// magenta picture"): a tint of +150 is a white 0.05 above the Planckian locus in CIE 1960 uv +/// (a greener light, which the white balance compensates towards magenta), −150 as far below it. +/// D65 reads about +10. +pub const TINT_PER_DUV: f64 = 3000.0; + +/// CIE 1960 uv of a Planckian radiator at `t` kelvin, 1 000–15 000 K (M. Krystek, "An algorithm +/// to calculate correlated colour temperature", Color Research & Application 10(1), 1985: a +/// rational approximation accurate to better than 1e-4 in uv). +pub fn planck_uv(t: f64) -> (f64, f64) { + let t = if t.is_finite() { t.clamp(1_000.0, 15_000.0) } else { 6_500.0 }; + let u = (0.860_117_757 + 1.541_182_54e-4 * t + 1.286_412_12e-7 * t * t) / (1.0 + 8.424_202_35e-4 * t + 7.081_451_63e-7 * t * t); + let v = (0.317_398_726 + 4.228_062_45e-5 * t + 4.204_816_91e-8 * t * t) / (1.0 - 2.897_418_16e-5 * t + 1.614_560_53e-7 * t * t); + (u, v) +} + +/// The unit normal of the Planckian locus at `t` in uv, pointing towards green (larger v). +fn planck_normal(t: f64) -> (f64, f64) { + let t = if t.is_finite() { t.clamp(1_000.0, 15_000.0) } else { 6_500.0 }; + let (a, b) = (planck_uv((t - 1.0).max(1_000.0)), planck_uv((t + 1.0).min(15_000.0))); + let (du, dv) = (b.0 - a.0, b.1 - a.1); + let len = (du * du + dv * dv).sqrt(); + if len < 1e-15 { + return (0.0, 1.0); + } + // perpendicular to the tangent; the locus runs right to left as t grows, so (dv, -du) points up + let n = (-dv / len, du / len); + if n.1 < 0.0 { (-n.0, -n.1) } else { n } +} + +fn uv_to_xy((u, v): (f64, f64)) -> Option<(f64, f64)> { + let d = 2.0 * u - 8.0 * v + 4.0; + (d.abs() > 1e-12).then(|| (3.0 * u / d, 2.0 * v / d)) +} + +fn xy_to_uv((x, y): (f64, f64)) -> Option<(f64, f64)> { + let d = -2.0 * x + 12.0 * y + 3.0; + (d.abs() > 1e-12).then(|| (4.0 * x / d, 6.0 * y / d)) +} + +/// The white chromaticity of a colour temperature (kelvin, 1 000–15 000 on the Planckian locus) +/// and tint ([`TINT_PER_DUV`] units; positive makes the developed picture more magenta, negative +/// greener, ±150 is the useful range). +pub fn temperature_tint_to_xy(temperature: f64, tint: f64) -> Option<(f64, f64)> { + if !temperature.is_finite() || !tint.is_finite() { + return None; + } + let t = temperature.clamp(1_000.0, 15_000.0); + let (u, v) = planck_uv(t); + let (nu, nv) = planck_normal(t); + let duv = tint.clamp(-1_000.0, 1_000.0) / TINT_PER_DUV; + uv_to_xy((u + nu * duv, v + nv * duv)) +} + +/// The colour temperature and tint of a white chromaticity (the inverse of +/// [`temperature_tint_to_xy`]): the nearest point of the Planckian locus in uv, searched in +/// reciprocal temperature. +pub fn xy_to_temperature_tint(xy: (f64, f64)) -> Option<(f64, f64)> { + let (u, v) = xy_to_uv(xy)?; + if !u.is_finite() || !v.is_finite() { + return None; + } + let dist = |mired: f64| { + let (pu, pv) = planck_uv(1e6 / mired); + (u - pu).powi(2) + (v - pv).powi(2) + }; + // coarse scan over 1 000–15 000 K, then golden-section refinement + let (lo_m, hi_m) = (1e6 / 15_000.0, 1e6 / 1_000.0); + let steps = 200; + let mut best = lo_m; + for i in 0..=steps { + let m = lo_m + (hi_m - lo_m) * f64::from(i) / f64::from(steps); + if dist(m) < dist(best) { + best = m; + } + } + let step = (hi_m - lo_m) / f64::from(steps); + let (mut a, mut b) = ((best - step).max(lo_m), (best + step).min(hi_m)); + let g = (5f64.sqrt() - 1.0) / 2.0; + for _ in 0..80 { + let (c, d) = (b - g * (b - a), a + g * (b - a)); + if dist(c) < dist(d) { b = d } else { a = c } + } + let t = 1e6 / ((a + b) / 2.0); + let (pu, pv) = planck_uv(t); + let (nu, nv) = planck_normal(t); + let duv = (u - pu) * nu + (v - pv) * nv; + Some((t, duv * TINT_PER_DUV)) +} + /// Colour temperature of an EXIF LightSource code used as a DNG CalibrationIlluminant. pub fn illuminant_kelvin(code: u64) -> Option { Some(match code { @@ -372,6 +462,32 @@ mod tests { } } + #[test] + fn temperature_and_tint_round_trip() { + // the Planckian locus at 6504 K is within 0.0035 of D65 in uv (D65 sits above it) + let xy = temperature_tint_to_xy(6504.0, 0.0).unwrap(); + assert!((xy.0 - D65.0).abs() < 0.005 && (xy.1 - D65.1).abs() < 0.01, "{xy:?}"); + // illuminant A is a Planckian radiator at 2856 K + let a = temperature_tint_to_xy(2856.0, 0.0).unwrap(); + assert!((a.0 - 0.44757).abs() < 5e-4 && (a.1 - 0.40745).abs() < 5e-4, "{a:?}"); + for t in [2000.0, 3200.0, 5600.0, 9000.0, 14_000.0] { + for tint in [-100.0, 0.0, 37.0, 150.0] { + let (t2, tint2) = xy_to_temperature_tint(temperature_tint_to_xy(t, tint).unwrap()).unwrap(); + assert!((t2 - t).abs() / t < 1e-4, "{t} {tint}: {t2}"); + assert!((tint2 - tint).abs() < 1e-3, "{t} {tint}: {tint2}"); + } + } + // positive tint is a greener white (above the locus), so the picture turns magenta + let m = temperature_tint_to_xy(5000.0, 50.0).unwrap(); + let g = temperature_tint_to_xy(5000.0, -50.0).unwrap(); + assert!(m.1 > g.1); + // D65 sits above the locus: about +10 + let (_, d65_tint) = xy_to_temperature_tint(D65).unwrap(); + assert!((8.0..12.0).contains(&d65_tint), "{d65_tint}"); + assert!(temperature_tint_to_xy(f64::NAN, 0.0).is_none()); + assert!(xy_to_temperature_tint((0.0, 0.0)).is_some_and(|(t, _)| t.is_finite())); + } + #[test] fn hostile_calibrations_are_errors() { let zero = Calibration { color_matrix: [Some([[0.0; 3]; 3]), None], as_shot_neutral: Some([1.0; 3]), ..Default::default() }; diff --git a/crates/dng/src/demosaic.rs b/crates/dng/src/demosaic.rs index 911236929b..e35909b7dc 100644 --- a/crates/dng/src/demosaic.rs +++ b/crates/dng/src/demosaic.rs @@ -13,6 +13,8 @@ #[cfg(feature = "threads")] use rayon::prelude::*; +use crate::options::DemosaicQuality; + /// A CFA pattern: `rows` × `cols` colour indices (0 red, 1 green, 2 blue), origin at the top left /// of the image. #[derive(Clone, Debug, PartialEq, Eq)] @@ -71,12 +73,17 @@ fn padded(plane: &[f32], w: usize, h: usize) -> (Vec, usize) { /// Demosaic `plane` (`w` × `h`, one sample per pixel) into interleaved RGB. pub fn demosaic(plane: &[f32], w: usize, h: usize, pattern: &Pattern) -> Vec { + demosaic_with(plane, w, h, pattern, DemosaicQuality::High) +} + +/// [`demosaic`] at a chosen quality: [`DemosaicQuality::Fast`] is bilinear for every pattern. +pub fn demosaic_with(plane: &[f32], w: usize, h: usize, pattern: &Pattern, quality: DemosaicQuality) -> Vec { let mut out = vec![0f32; w * h * 3]; if w == 0 || h == 0 || plane.len() < w * h { return out; } let (p, pw) = padded(plane, w, h); - let bayer = pattern.is_bayer() && w >= 3 && h >= 3; + let bayer = quality == DemosaicQuality::High && pattern.is_bayer() && w >= 3 && h >= 3; let row = |y: usize, o: &mut [f32]| { if bayer { malvar_row(&p, pw, y, w, pattern, o); @@ -91,6 +98,126 @@ pub fn demosaic(plane: &[f32], w: usize, h: usize, pattern: &Pattern) -> Vec { + rows: [&'a [f32]; 5], + w: usize, + /// Colour of even and odd columns of this row. + colors: [u8; 2], + high: bool, +} + +impl<'a> BayerRow<'a> { + /// Row `y` of the `w` × `h` plane (`w`, `h` ≥ 3 and a 2 × 2 Bayer `pat`). + pub(crate) fn new(plane: &'a [f32], w: usize, h: usize, pat: &Pattern, quality: DemosaicQuality, y: usize) -> BayerRow<'a> { + let row = |dy: isize| { + let sy = mirror(y as isize + dy, h); + plane.get(sy * w..sy * w + w).unwrap_or(&[]) + }; + BayerRow { rows: [row(-2), row(-1), row(0), row(1), row(2)], w, colors: [pat.color(0, y), pat.color(1, y)], high: quality == DemosaicQuality::High } + } + + /// The RGB of columns `x0..x1`, in order, passed to `f`: the interior through 5 × 5 windows + /// sliding along the five rows, the two columns at either edge through [`BayerRow::pixel`]. + #[inline] + pub(crate) fn span(&mut self, x0: usize, x1: usize, mut f: impl FnMut([f32; 3])) { + let (lo, hi) = (2.min(self.w), self.w.saturating_sub(2).max(2)); + let mut rgb = [0f32; 3]; + let mut x = x0; + while x < x1.min(lo) { + self.pixel(x, &mut rgb); + f(rgb); + x += 1; + } + let end = x1.min(hi); + if x < end { + let [r0, r1, r2, r3, r4] = self.rows; + let win = |r: &'a [f32]| r.get(x - 2..end + 2).unwrap_or(&[]).windows(5); + let rows = win(r0).zip(win(r1)).zip(win(r2)).zip(win(r3)).zip(win(r4)); + let (high, colors) = (self.high, self.colors); + for ((((a, b), c), d), e) in rows { + let at = |dx: isize, dy: isize| { + let r = match dy { + -2 => a, + -1 => b, + 0 => c, + 1 => d, + _ => e, + }; + r.get((dx + 2) as usize).copied().unwrap_or(0.0) + }; + let me = colors[x & 1]; + let right_is_red = colors[(x + 1) & 1] == 0; + f(if high { malvar_px(at, me, right_is_red) } else { bilinear_px(at, me, right_is_red) }); + x += 1; + } + } + while x < x1 { + self.pixel(x, &mut rgb); + f(rgb); + x += 1; + } + } + + /// The RGB of column `x`. + #[inline] + pub(crate) fn pixel(&mut self, x: usize, out: &mut [f32; 3]) { + let me = self.colors[x & 1]; + let right_is_red = self.colors[(x + 1) & 1] == 0; + let rows = &self.rows; + if x >= 2 && x + 2 < self.w { + // the 5 × 5 window, bounds checked once per row + let win = rows.map(|r| r.get(x - 2..x + 3).and_then(|s| <[f32; 5]>::try_from(s).ok()).unwrap_or([0.0; 5])); + let at = |dx: isize, dy: isize| win[(dy + 2) as usize][(dx + 2) as usize]; + *out = if self.high { malvar_px(at, me, right_is_red) } else { bilinear_px(at, me, right_is_red) }; + } else { + let w = self.w; + let at = |dx: isize, dy: isize| rows[(dy + 2) as usize].get(mirror(x as isize + dx, w)).copied().unwrap_or(0.0); + *out = if self.high { malvar_px(at, me, right_is_red) } else { bilinear_px(at, me, right_is_red) }; + } + } +} + +/// Malvar–He–Cutler at one pixel of colour `me` (the expressions of [`malvar_row`]). +#[inline(always)] +fn malvar_px(at: impl Fn(isize, isize) -> f32, me: u8, right_is_red: bool) -> [f32; 3] { + let c = at(0, 0); + let diag = at(-1, -1) + at(1, -1) + at(-1, 1) + at(1, 1); + if me == 1 { + let h1 = at(-1, 0) + at(1, 0); + let h2 = at(-2, 0) + at(2, 0); + let v1 = at(0, -1) + at(0, 1); + let v2 = at(0, -2) + at(0, 2); + let along_row = (5.0 * c + 4.0 * h1 - h2 - diag + 0.5 * v2) / 8.0; + let along_col = (5.0 * c + 4.0 * v1 - v2 - diag + 0.5 * h2) / 8.0; + if right_is_red { [along_row, c, along_col] } else { [along_col, c, along_row] } + } else { + let cross1 = at(-1, 0) + at(1, 0) + at(0, -1) + at(0, 1); + let cross2 = at(-2, 0) + at(2, 0) + at(0, -2) + at(0, 2); + let g = (4.0 * c + 2.0 * cross1 - cross2) / 8.0; + let other = (6.0 * c + 2.0 * diag - 1.5 * cross2) / 8.0; + if me == 0 { [c, g, other] } else { [other, g, c] } + } +} + +/// Bilinear at one Bayer pixel: the sums of [`bilinear_row`]'s 3 × 3 window in its order. +#[inline(always)] +fn bilinear_px(at: impl Fn(isize, isize) -> f32, me: u8, right_is_red: bool) -> [f32; 3] { + let c = at(0, 0); + if me == 1 { + let along_row = (0.0 + at(-1, 0) + at(1, 0)) / 2.0; + let along_col = (0.0 + at(0, -1) + at(0, 1)) / 2.0; + if right_is_red { [along_row, c, along_col] } else { [along_col, c, along_row] } + } else { + let g = (0.0 + at(0, -1) + at(-1, 0) + at(1, 0) + at(0, 1)) / 4.0; + let other = (0.0 + at(-1, -1) + at(1, -1) + at(-1, 1) + at(1, 1)) / 4.0; + if me == 0 { [c, g, other] } else { [other, g, c] } + } +} + fn malvar_row(p: &[f32], pw: usize, y: usize, w: usize, pat: &Pattern, out: &mut [f32]) { let py = y + PAD; for x in 0..w { diff --git a/crates/dng/src/lib.rs b/crates/dng/src/lib.rs index c0b3e9e0ba..e136d1d321 100644 --- a/crates/dng/src/lib.rs +++ b/crates/dng/src/lib.rs @@ -17,14 +17,19 @@ pub mod color; pub mod demosaic; pub mod ljpeg; +pub mod opcode; +pub mod options; pub mod process; pub mod raw; +#[cfg(any(test, feature = "writer"))] +pub mod reference; pub mod tiff; #[cfg(any(test, feature = "writer"))] pub mod writer; use color::{Calibration, M3, V3}; use demosaic::Pattern; +pub use options::{DemosaicQuality, Highlights, OutputSpace, ProcessOptions, WhiteBalance}; use process::Params; use raw::{PHOTOMETRIC_CFA, RawImage}; use tiff::{Entry, Ifd, Tiff, tags}; @@ -89,7 +94,7 @@ pub struct Info { pub samples_per_pixel: usize, /// Colour filter array data (else LinearRaw). pub cfa: bool, - /// "uncompressed" or "lossless JPEG". + /// "uncompressed", "lossless JPEG" or "lossy JPEG". pub compression: &'static str, pub dng_version: [u8; 4], /// CinemaDNG FrameRate (numerator, denominator). @@ -213,6 +218,11 @@ impl<'a> Dng<'a> { /// The processing tags of the raw image. pub fn params(&self) -> Result { + self.params_with(true) + } + + /// [`Dng::params`], reading the opcode lists only when `opcodes` is set. + fn params_with(&self, opcodes: bool) -> Result { let uint_pair = |tag| -> Result> { match self.raw_tag(tag) { Some(e) => match self.tiff.uints(e)?[..] { @@ -260,9 +270,21 @@ impl<'a> Dng<'a> { crop_origin: pair(tags::DEFAULT_CROP_ORIGIN)?, crop_size: pair(tags::DEFAULT_CROP_SIZE)?, cfa: self.cfa()?, + opcodes: opcode::Lists { + list1: if opcodes { self.opcode_list(tags::OPCODE_LIST_1)? } else { Vec::new() }, + list2: if opcodes { self.opcode_list(tags::OPCODE_LIST_2)? } else { Vec::new() }, + list3: if opcodes { self.opcode_list(tags::OPCODE_LIST_3)? } else { Vec::new() }, + }, }) } + fn opcode_list(&self, tag: u16) -> Result> { + match self.raw_ifd().get(tag) { + Some(e) => opcode::parse(self.tiff.value_bytes(e)?), + None => Ok(Vec::new()), + } + } + /// The colour tags. pub fn calibration(&self) -> Result { let illum = |tag| self.main_tag(tag).map(|e| self.tiff.uint(e)).transpose(); @@ -285,7 +307,7 @@ impl<'a> Dng<'a> { /// Size, depth and CinemaDNG metadata. pub fn info(&self) -> Result { let l = raw::layout(&self.tiff, self.raw_ifd())?; - let p = self.params()?; + let p = self.params_with(false)?; let (_, _, aw, ah) = process::active_rect(&p, l.width, l.height)?; let (_, _, width, height) = process::crop_rect(&p, aw, ah); let dng_version = self @@ -315,7 +337,11 @@ impl<'a> Dng<'a> { bits: l.bits, samples_per_pixel: l.samples_per_pixel, cfa: p.cfa.is_some(), - compression: if l.compression == 7 { "lossless JPEG" } else { "uncompressed" }, + compression: match l.compression { + 7 => "lossless JPEG", + 34892 => "lossy JPEG", + _ => "uncompressed", + }, dng_version, frame_rate, timecode, @@ -329,12 +355,39 @@ impl<'a> Dng<'a> { raw::read(&self.tiff, self.raw_ifd()) } - /// Develop the main image to scene-linear Rec. 709. + /// The as-shot white balance as (temperature in kelvin, tint), for a settings panel to start + /// from ([`WhiteBalance::Custom`] with these values develops like [`WhiteBalance::AsShot`]). + /// `None` for monochrome files and files without a ColorMatrix. + pub fn as_shot_temperature_tint(&self) -> Result> { + let calibration = self.calibration()?; + if calibration.color_matrix.iter().all(Option::is_none) || self.cfa()?.is_none() && raw::layout(&self.tiff, self.raw_ifd())?.samples_per_pixel == 1 { + return Ok(None); + } + Ok(color::xy_to_temperature_tint(calibration.solve()?.white_xy)) + } + + /// Develop the main image to scene-linear Rec. 709 as the file says + /// ([`ProcessOptions::default`]). pub fn develop(&self) -> Result { + self.develop_with(&ProcessOptions::default()) + } + + /// Develop the main image with raw settings. + pub fn develop_with(&self, o: &ProcessOptions) -> Result { let info = self.info()?; - let p = self.params()?; - let raw = self.raw()?; - let calibration = self.calibration()?; + let p = self.params_with(o.apply_opcodes)?; + let mut raw = self.raw()?; + let preview = process::preview(&p, o, &info); + process::apply_list1(&mut raw, &p, preview)?; + let mut calibration = self.calibration()?; + calibration.baseline_exposure += o.exposure_stops(); + if let WhiteBalance::Custom { temperature, tint } = o.white_balance + && calibration.color_matrix.iter().any(Option::is_some) + { + let xy = color::temperature_tint_to_xy(temperature, tint).ok_or(Error::Unsupported("a non-finite white balance".into()))?; + calibration.as_shot_neutral = None; + calibration.as_shot_white_xy = Some(xy); + } let solved = if p.cfa.is_none() && raw.samples_per_pixel == 1 { // monochrome: no colour, only the exposure let e = 2f64.powf(calibration.baseline_exposure.clamp(-16.0, 16.0)); @@ -342,7 +395,7 @@ impl<'a> Dng<'a> { } else { calibration.solve()? }; - let d = process::develop(&raw, &p, &solved.neutral, &solved.matrix)?; + let d = process::develop_with(&raw, &p, &solved.neutral, &solved.matrix, o)?; Ok(Image { info, width: d.width, height: d.height, rgba: d.rgba }) } } @@ -352,6 +405,11 @@ pub fn decode(bytes: &[u8]) -> Result { Dng::parse(bytes)?.develop() } +/// Parse and develop a DNG file with raw settings. +pub fn decode_with(bytes: &[u8], o: &ProcessOptions) -> Result { + Dng::parse(bytes)?.develop_with(o) +} + /// Parse a DNG file's metadata. pub fn probe(bytes: &[u8]) -> Result { Dng::parse(bytes)?.info() diff --git a/crates/dng/src/ljpeg.rs b/crates/dng/src/ljpeg.rs index 15e39fa827..25cf4da8fd 100644 --- a/crates/dng/src/ljpeg.rs +++ b/crates/dng/src/ljpeg.rs @@ -27,6 +27,10 @@ pub struct Decoded { } const LOOKUP_BITS: u32 = 9; +/// Bits of the combined code + difference lookup. +const FULL_BITS: u32 = 13; +/// A combined lookup entry: bit 8 set, the bits used in 0–7, the difference in 16–31. +const FULL_VALID: u32 = 0x100; /// A Huffman table for decoding (T.81 Annex C and F.2.2.3). #[derive(Clone)] @@ -38,6 +42,10 @@ struct Huffman { valptr: [i32; 17], mincode: [i32; 17], values: Vec, + /// For a code and its extra bits that fit in FULL_BITS together (difference categories + /// 0–15): `FULL_VALID | used bits | difference << 16`, indexed by the next FULL_BITS bits; + /// 0 means "decode the slow way". + full: Vec, } impl Huffman { @@ -96,7 +104,28 @@ impl Huffman { } } } - Ok(Huffman { fast, maxcode, valptr, mincode, values }) + let mut full = vec![0u32; 1 << FULL_BITS]; + for (&c, (&s, &ssss)) in codes.iter().zip(sizes.iter().zip(&values)) { + let (s, ssss) = (u32::from(s), u32::from(ssss)); + if ssss >= 16 || s + ssss > FULL_BITS { + continue; + } + let total = s + ssss; + for v in 0..1u32 << ssss { + let diff = match ssss { + 0 => 0, + k if v < 1 << (k - 1) => v as i32 - (1 << k) + 1, + _ => v as i32, + }; + let entry = FULL_VALID | total | (u32::from(diff as i16 as u16) << 16); + let shift = FULL_BITS - total; + let base = (((c << ssss) | v) << shift) as usize; + for f in full.iter_mut().skip(base).take(1 << shift) { + *f = entry; + } + } + } + Ok(Huffman { fast, maxcode, valptr, mincode, values, full }) } } @@ -116,6 +145,23 @@ impl<'a> Bits<'a> { } fn fill(&mut self) { + // fast path: eight bytes without a 0xFF (no stuffing, no marker) go in at once + if self.marker.is_none() + && let Some(b) = self.data.get(self.pos..self.pos + 8) + { + let chunk = u64::from_be_bytes([b[0], b[1], b[2], b[3], b[4], b[5], b[6], b[7]]); + let inv = !chunk; + let has_ff = inv.wrapping_sub(0x0101_0101_0101_0101) & !inv & 0x8080_8080_8080_8080 != 0; + if !has_ff { + let take = (64 - self.n) / 8; + if take > 0 { + self.acc |= (chunk >> (64 - 8 * take)) << (64 - self.n - 8 * take); + self.pos += take as usize; + self.n += 8 * take; + return; + } + } + } while self.n <= 56 { let mut byte = 0u8; if self.marker.is_none() { @@ -342,10 +388,16 @@ struct ScanParams { restart: usize, } -/// Read one difference value: a Huffman-coded category followed by its extra bits. -fn read_diff(bits: &mut Bits, h: &Huffman) -> Result { - let ssss = u32::from(bits.decode(h)?); - Ok(match ssss { +/// The next difference (F.2.2.1 and H.1.2.2): category and extra bits from the combined table +/// when they fit, else symbol by symbol. +#[inline(always)] +fn difference(bits: &mut Bits, h: &Huffman) -> Result { + let e = h.full.get(bits.peek(FULL_BITS) as usize).copied().unwrap_or(0); + if e & FULL_VALID != 0 { + bits.skip(e & 0xFF); + return Ok(i32::from((e >> 16) as u16 as i16)); + } + Ok(match u32::from(bits.decode(h)?) { 0 => 0, 16 => 32768, s => { @@ -373,26 +425,17 @@ fn scan(data: &[u8], start: usize, hufs: &[&Huffman], params: &ScanParams) -> Re bits.restart()?; line_in_interval = 0; } - let first_line = line_in_interval == 0; - let (above, cur) = if y == 0 { - (None, &mut out[..row_len]) - } else { - let (a, c) = out.split_at_mut(y * row_len); - (Some(&a[(y - 1) * row_len..]), &mut c[..row_len]) - }; - for x in 0..width { - for (c, h) in hufs.iter().enumerate() { - let i = x * nc + c; - let diff = read_diff(&mut bits, h)?; - let px = match above { - _ if first_line && x == 0 => initial, - _ if first_line => i32::from(cur[i - nc]), - Some(a) if x == 0 => i32::from(a[i]), - Some(a) => predict(predictor, i32::from(cur[i - nc]), i32::from(a[i]), i32::from(a[i - nc])), - None => initial, - }; - cur[i] = ((px + diff) & 0xFFFF) as u16; - } + let (above, cur) = out.split_at_mut(y * row_len); + let cur = cur.get_mut(..row_len).ok_or(Error::TooLarge)?; + let above = if line_in_interval == 0 { None } else { above.get((y - 1) * row_len..) }; + match predictor { + 1 => line::<1>(&mut bits, hufs, cur, above, initial)?, + 2 => line::<2>(&mut bits, hufs, cur, above, initial)?, + 3 => line::<3>(&mut bits, hufs, cur, above, initial)?, + 4 => line::<4>(&mut bits, hufs, cur, above, initial)?, + 5 => line::<5>(&mut bits, hufs, cur, above, initial)?, + 6 => line::<6>(&mut bits, hufs, cur, above, initial)?, + _ => line::<7>(&mut bits, hufs, cur, above, initial)?, } line_in_interval += 1; } @@ -404,7 +447,52 @@ fn scan(data: &[u8], start: usize, hufs: &[&Huffman], params: &ScanParams) -> Re Ok(out) } +/// One line of samples (H.1.2.1): the first line of a restart interval (`above` is `None`) +/// predicts from the left with `initial` for its first sample; other lines predict their first +/// sample from above and the rest with predictor `P`. +#[inline(always)] +fn line(bits: &mut Bits, hufs: &[&Huffman], cur: &mut [u16], above: Option<&[u16]>, initial: i32) -> Result<()> { + let nc = hufs.len(); + let row_len = cur.len(); + let above = match above { + Some(a) => Some(a.get(..row_len).ok_or(Error::TooLarge)?), + None => None, + }; + for (c, h) in hufs.iter().enumerate() { + let d = difference(bits, h)?; + let px = match above { + Some(a) => i32::from(a[c]), + None => initial, + }; + cur[c] = ((px + d) & 0xFFFF) as u16; + } + let mut i = nc; + match above { + None => { + while i < row_len { + for h in hufs { + let d = difference(bits, h)?; + cur[i] = ((i32::from(cur[i - nc]) + d) & 0xFFFF) as u16; + i += 1; + } + } + } + Some(a) => { + while i < row_len { + for h in hufs { + let d = difference(bits, h)?; + let px = predict(P, i32::from(cur[i - nc]), i32::from(a[i]), i32::from(a[i - nc])); + cur[i] = ((px + d) & 0xFFFF) as u16; + i += 1; + } + } + } + } + Ok(()) +} + /// Table H.1 predictors (`ra` left, `rb` above, `rc` above left). +#[inline(always)] pub fn predict(p: u8, ra: i32, rb: i32, rc: i32) -> i32 { match p { 1 => ra, diff --git a/crates/dng/src/opcode.rs b/crates/dng/src/opcode.rs new file mode 100644 index 0000000000..03bb86a55b --- /dev/null +++ b/crates/dng/src/opcode.rs @@ -0,0 +1,886 @@ +//! DNG opcode lists (DNG 1.7.1.0 chapter 7, "Opcode List Processing"): OpcodeList1 on the raw +//! samples as stored, OpcodeList2 on the linear reference values (after linearisation, black and +//! white levels, before white balance), OpcodeList3 on the demosaiced picture. +//! +//! An opcode list is big-endian whatever the file's byte order: a count, then per opcode its ID, +//! the DNG version that defined it, flags (bit 0 optional, bit 1 skippable for preview quality) +//! and the size of its parameters. Implemented: WarpRectilinear (1), WarpFisheye (2), +//! FixVignetteRadial (3), FixBadPixelsConstant (4), FixBadPixelsList (5), TrimBounds (6), +//! MapTable (7), MapPolynomial (8), GainMap (9), DeltaPerRow / Column (10, 11), ScalePerRow / +//! Column (12, 13) and WarpRectilinear2 (14, with the skip rule for a following warp). An unknown +//! opcode, or one from a DNG version newer than 1.7.1.0, is skipped when it is optional and an +//! error when it is not. Values are clipped to the stage's logical range after every opcode +//! (0–65535 for OpcodeList1, 0–1 after). +//! +//! Interpretations where the specification leaves room (recorded in the crate README): +//! coordinates of OpcodeList1 are those of the stored raw image, of OpcodeList2 and 3 those of the +//! ActiveArea; GainMap positions are relative to the whole image of the stage (row / height, +//! column / width); Delta and MapPolynomial values are in the 0–1 range at every stage (scaled by +//! 65535 on OpcodeList1's integer samples), MapTable indexes 16-bit values (0–1 scaled by 65535 +//! after OpcodeList1); warps resample with a Catmull–Rom cubic and are refused on CFA mosaics, +//! as is TrimBounds before OpcodeList3. +//! +//! Every count and size is checked against the bytes of the parameter area before anything is +//! allocated, and per-opcode work is bounded by the image size. + +use crate::{Error, Result}; + +/// The newest DNG version whose opcodes this reader knows (1.7.1.0). +const KNOWN_VERSION: u32 = 0x0107_0100; + +/// Most opcodes in one list (each may touch the whole picture). +pub const MAX_OPCODES: usize = 256; + +const OPTIONAL: u32 = 1; +const PREVIEW_SKIPPABLE: u32 = 2; + +/// One opcode of a list. +#[derive(Clone, Debug, PartialEq)] +pub struct Opcode { + pub id: u32, + /// The DNG version that defined it, as four bytes in a big-endian u32 (1.3.0.0 = 0x01030000). + pub version: u32, + pub flags: u32, + pub params: Vec, +} + +impl Opcode { + pub fn optional(&self) -> bool { + self.flags & OPTIONAL != 0 + } + + pub fn preview_skippable(&self) -> bool { + self.flags & PREVIEW_SKIPPABLE != 0 + } +} + +/// The opcode lists of a raw IFD. +#[derive(Clone, Debug, Default, PartialEq)] +pub struct Lists { + pub list1: Vec, + pub list2: Vec, + pub list3: Vec, +} + +/// Parse an opcode list (big-endian). +pub fn parse(b: &[u8]) -> Result> { + let bad = || Error::Malformed("opcode list"); + let mut r = Reader::new(b); + let count = r.u32().map_err(|_| bad())? as usize; + // every opcode takes at least 16 bytes + if count > b.len() / 16 || count > MAX_OPCODES { + return Err(bad()); + } + let mut out = Vec::with_capacity(count); + for _ in 0..count { + let id = r.u32().map_err(|_| bad())?; + let version = r.u32().map_err(|_| bad())?; + let flags = r.u32().map_err(|_| bad())?; + let len = r.u32().map_err(|_| bad())? as usize; + let params = r.bytes(len).map_err(|_| bad())?.to_vec(); + out.push(Opcode { id, version, flags, params }); + } + Ok(out) +} + +/// Big-endian parameter reader. +struct Reader<'a> { + b: &'a [u8], + pos: usize, +} + +impl<'a> Reader<'a> { + fn new(b: &'a [u8]) -> Self { + Reader { b, pos: 0 } + } + + fn bytes(&mut self, n: usize) -> Result<&'a [u8]> { + let end = self.pos.checked_add(n).ok_or(Error::Malformed("opcode parameters"))?; + let s = self.b.get(self.pos..end).ok_or(Error::Malformed("opcode parameters too short"))?; + self.pos = end; + Ok(s) + } + + fn u32(&mut self) -> Result { + let s = self.bytes(4)?; + Ok(u32::from_be_bytes([s[0], s[1], s[2], s[3]])) + } + + fn f32(&mut self) -> Result { + let v = f32::from_bits(self.u32()?); + if v.is_finite() { Ok(v) } else { Err(Error::Malformed("non-finite opcode parameter")) } + } + + fn f64(&mut self) -> Result { + let s = self.bytes(8)?; + let v = f64::from_be_bytes([s[0], s[1], s[2], s[3], s[4], s[5], s[6], s[7]]); + if v.is_finite() { Ok(v) } else { Err(Error::Malformed("non-finite opcode parameter")) } + } + + fn remaining(&self) -> usize { + self.b.len().saturating_sub(self.pos) + } + + /// A count of items of `size` bytes each that must fit in what is left. + fn count(&mut self, size: usize) -> Result { + let n = self.u32()? as usize; + if n.saturating_mul(size) > self.remaining() { + return Err(Error::Malformed("opcode count larger than its parameters")); + } + Ok(n) + } +} + +/// The picture an opcode list works on: `planes` interleaved samples per pixel. +pub struct Image<'a> { + pub data: &'a mut Vec, + pub width: usize, + pub height: usize, + pub planes: usize, + /// The stage's logical maximum: 65535 for OpcodeList1 (integer samples, rounded after every + /// opcode), 1 for OpcodeList2 and 3. + pub max: f32, + /// A CFA mosaic (one sample per pixel before the demosaic): warps are refused. + pub mosaic: bool, + /// TrimBounds is allowed (OpcodeList3 only: the picture may change size). + pub allow_trim: bool, + /// Origin of the picture after TrimBounds, relative to where it started. + pub origin: (usize, usize), +} + +impl Image<'_> { + fn integer(&self) -> bool { + self.max > 1.0 + } + + fn clip(&mut self) { + let (max, int) = (self.max, self.integer()); + for v in self.data.iter_mut() { + let c = v.clamp(0.0, max); + *v = if int { + c.round() + } else if c.is_nan() { + 0.0 + } else { + c + }; + } + } +} + +/// Apply `list`. `preview` skips the opcodes marked as only needed for full quality. +pub fn apply(list: &[Opcode], img: &mut Image, preview: bool) -> Result<()> { + let mut skip_warp = false; + for op in list { + let skip_this = std::mem::take(&mut skip_warp) && matches!(op.id, 1 | 2); + if skip_this || (preview && op.preview_skippable()) { + continue; + } + let known = op.version <= KNOWN_VERSION && (1..=14).contains(&op.id); + if !known { + if op.optional() { + continue; + } + return Err(Error::Unsupported(format!("DNG opcode {} (version {:08x})", op.id, op.version))); + } + match run(op, img) { + Ok(()) => {} + // an optional opcode this reader cannot apply here is left out + Err(Error::Unsupported(_)) if op.optional() => continue, + Err(e) => return Err(e), + } + img.clip(); + if op.id == 14 && op.optional() { + skip_warp = true; + } + } + Ok(()) +} + +fn run(op: &Opcode, img: &mut Image) -> Result<()> { + let mut r = Reader::new(&op.params); + match op.id { + 1 | 2 | 14 => warp(op.id, &mut r, img), + 3 => vignette(&mut r, img), + 4 => { + let constant = r.u32()?; + let phase = r.u32()?; + let mut bad = Vec::new(); + let c = if img.integer() { constant as f32 } else { constant as f32 / 65_535.0 }; + for (i, px) in img.data.iter().enumerate() { + if *px == c { + bad.push(i); + } + } + fix_bad(img, phase, &bad) + } + 5 => { + let phase = r.u32()?; + let points = r.u32()? as usize; + let rects = r.u32()? as usize; + if points.saturating_mul(8).saturating_add(rects.saturating_mul(16)) > r.remaining() { + return Err(Error::Malformed("FixBadPixelsList counts larger than its parameters")); + } + let (w, h) = (img.width, img.height); + let mut mask = vec![false; w * h]; + for _ in 0..points { + let (y, x) = (r.u32()? as usize, r.u32()? as usize); + if x < w && y < h { + mask[y * w + x] = true; + } + } + for _ in 0..rects { + let (t, l, b, rt) = (r.u32()? as usize, r.u32()? as usize, r.u32()? as usize, r.u32()? as usize); + for y in t.min(h)..b.min(h) { + for x in l.min(w)..rt.min(w) { + mask[y * w + x] = true; + } + } + } + let bad: Vec = mask.iter().enumerate().filter_map(|(i, &b)| b.then_some(i)).collect(); + fix_bad(img, phase, &bad) + } + 6 => trim(&mut r, img), + _ => area_op(op.id, &mut r, img), + } +} + +/// The area fields of opcodes 7–13, clamped to the picture. +struct Area { + top: usize, + left: usize, + bottom: usize, + right: usize, + plane: usize, + planes: usize, + row_pitch: usize, + col_pitch: usize, +} + +impl Area { + fn read(r: &mut Reader, img: &Image) -> Result { + let mut v = [0usize; 8]; + for x in &mut v { + *x = r.u32()? as usize; + } + let [top, left, bottom, right, plane, planes, row_pitch, col_pitch] = v; + if row_pitch == 0 || col_pitch == 0 { + return Err(Error::Malformed("opcode area with a zero pitch")); + } + let plane = plane.min(img.planes); + Ok(Area { + top: top.min(img.height), + left: left.min(img.width), + bottom: bottom.min(img.height), + right: right.min(img.width), + plane, + planes: planes.min(img.planes - plane), + row_pitch, + col_pitch, + }) + } + + /// Rows of the area, in order. + fn rows(&self) -> impl Iterator + use<> { + (self.top..self.bottom.max(self.top)).step_by(self.row_pitch) + } + + fn cols(&self) -> impl Iterator + use<> { + (self.left..self.right.max(self.left)).step_by(self.col_pitch) + } + + /// Visit every affected sample: f(row index in area, col index in area, x, y, plane index in + /// area, value). + fn each(&self, img: &mut Image, mut f: impl FnMut(usize, usize, usize, usize, usize, &mut f32)) { + let (w, np) = (img.width, img.planes); + for (ri, y) in self.rows().enumerate() { + for (ci, x) in self.cols().enumerate() { + for p in 0..self.planes { + if let Some(v) = img.data.get_mut((y * w + x) * np + self.plane + p) { + f(ri, ci, x, y, p, v); + } + } + } + } + } +} + +fn area_op(id: u32, r: &mut Reader, img: &mut Image) -> Result<()> { + let a = Area::read(r, img)?; + let scale = img.max; + match id { + 7 => { + let n = r.count(2)?; + if n == 0 { + return Err(Error::Malformed("MapTable without entries")); + } + let table: Vec = (0..n).map(|_| r.bytes(2).map(|b| f32::from(u16::from_be_bytes([b[0], b[1]])))).collect::>()?; + let last = table.last().copied().unwrap_or(0.0); + a.each(img, |_, _, _, _, _, v| { + let code = (*v / scale * 65_535.0).round().clamp(0.0, 65_535.0) as usize; + *v = table.get(code).copied().unwrap_or(last) / 65_535.0 * scale; + }); + } + 8 => { + let degree = r.u32()? as usize; + if degree > 8 { + return Err(Error::Malformed("MapPolynomial degree above 8")); + } + let c: Vec = (0..=degree).map(|_| r.f64()).collect::>()?; + a.each(img, |_, _, _, _, _, v| { + let x = f64::from(*v / scale); + let y = c.iter().rev().fold(0.0, |acc, k| acc * x + k); + *v = (y * f64::from(scale)) as f32; + }); + } + 9 => gain_map(&a, r, img)?, + 10..=13 => { + let n = r.count(4)?; + let vals: Vec = (0..n).map(|_| r.f32()).collect::>()?; + let per_row = matches!(id, 10 | 12); + let need = if per_row { a.rows().count() } else { a.cols().count() }; + if vals.len() < need { + return Err(Error::Malformed("opcode with fewer values than rows or columns")); + } + let delta = id <= 11; + a.each(img, |ri, ci, _, _, _, v| { + let k = vals.get(if per_row { ri } else { ci }).copied().unwrap_or(if delta { 0.0 } else { 1.0 }); + *v = if delta { *v + k * scale } else { *v * k }; + }); + } + _ => return Err(Error::Unsupported(format!("DNG opcode {id}"))), + } + Ok(()) +} + +fn gain_map(a: &Area, r: &mut Reader, img: &mut Image) -> Result<()> { + let points_v = r.u32()? as usize; + let points_h = r.u32()? as usize; + let spacing_v = r.f64()?; + let spacing_h = r.f64()?; + let origin_v = r.f64()?; + let origin_h = r.f64()?; + let map_planes = r.u32()? as usize; + let n = points_v.saturating_mul(points_h).saturating_mul(map_planes); + if points_v == 0 || points_h == 0 || map_planes == 0 || n.saturating_mul(4) > r.remaining() { + return Err(Error::Malformed("GainMap size")); + } + if (points_v > 1 && spacing_v <= 0.0) || (points_h > 1 && spacing_h <= 0.0) { + return Err(Error::Malformed("GainMap spacing")); + } + let gains: Vec = (0..n).map(|_| r.f32()).collect::>()?; + let (w, h) = (img.width.max(1) as f64, img.height.max(1) as f64); + let pos = |rel: f64, origin: f64, spacing: f64, points: usize| { + if points <= 1 { + return 0.0; + } + ((rel - origin) / spacing).clamp(0.0, (points - 1) as f64) + }; + let gain = |p: usize, mv: f64, mh: f64| -> f32 { + let mp = p.min(map_planes - 1); + let (v0, h0) = (mv.floor() as usize, mh.floor() as usize); + let (v1, h1) = ((v0 + 1).min(points_v - 1), (h0 + 1).min(points_h - 1)); + let (fv, fh) = (mv - v0 as f64, mh - h0 as f64); + let g = |v: usize, hh: usize| f64::from(gains.get((v * points_h + hh) * map_planes + mp).copied().unwrap_or(1.0)); + let top = g(v0, h0) * (1.0 - fh) + g(v0, h1) * fh; + let bottom = g(v1, h0) * (1.0 - fh) + g(v1, h1) * fh; + (top * (1.0 - fv) + bottom * fv) as f32 + }; + a.each(img, |_, _, x, y, p, v| { + let mv = pos(y as f64 / h, origin_v, spacing_v, points_v); + let mh = pos(x as f64 / w, origin_h, spacing_h, points_h); + *v *= gain(p, mv, mh); + }); + Ok(()) +} + +/// The optical centre in pixels and the distance to the farthest corner (chapter 7). +fn centre(img: &Image, cx: f64, cy: f64) -> (f64, f64, f64) { + let (x1, y1) = (img.width.saturating_sub(1) as f64, img.height.saturating_sub(1) as f64); + let (cx, cy) = (cx * x1, cy * y1); + let mx = cx.abs().max((x1 - cx).abs()); + let my = cy.abs().max((y1 - cy).abs()); + (cx, cy, (mx * mx + my * my).sqrt()) +} + +fn vignette(r: &mut Reader, img: &mut Image) -> Result<()> { + let k: Vec = (0..5).map(|_| r.f64()).collect::>()?; + let (cx, cy) = (r.f64()?, r.f64()?); + let (cx, cy, m) = centre(img, cx, cy); + if m <= 0.0 { + return Ok(()); + } + let (w, np) = (img.width, img.planes); + for (i, px) in img.data.chunks_exact_mut(np.max(1)).enumerate() { + let (x, y) = ((i % w) as f64, (i / w) as f64); + let r2 = ((x - cx).powi(2) + (y - cy).powi(2)) / (m * m); + let g = 1.0 + r2 * (k[0] + r2 * (k[1] + r2 * (k[2] + r2 * (k[3] + r2 * k[4])))); + for v in px { + *v = (f64::from(*v) * g) as f32; + } + } + Ok(()) +} + +/// Interpolate over `bad` pixels of a Bayer mosaic: the mean of the nearest good pixels of the +/// same colour (within 2, else 4, … up to 8 pixels). +fn fix_bad(img: &mut Image, phase: u32, bad: &[usize]) -> Result<()> { + if bad.is_empty() { + return Ok(()); + } + // a real sensor has few bad pixels; a list covering much of the picture is not one + if bad.len() > img.data.len() / 16 + 64 { + return Err(Error::Unsupported("bad pixel opcode covering much of the picture".into())); + } + if img.planes != 1 { + return Err(Error::Unsupported("bad pixel opcodes on a picture that is not a mosaic".into())); + } + // BayerPhase: the colour of the top-left pixel (0 red, 1 green in a red row, 2 green in a + // blue row, 3 blue) + let quad: [u8; 4] = match phase { + 0 => [0, 1, 1, 2], + 1 => [1, 0, 2, 1], + 2 => [1, 2, 0, 1], + 3 => [2, 1, 1, 0], + _ => return Err(Error::Malformed("BayerPhase")), + }; + let (w, h) = (img.width, img.height); + let mut is_bad = vec![false; w * h]; + for &i in bad { + if let Some(b) = is_bad.get_mut(i) { + *b = true; + } + } + let color = |x: usize, y: usize| quad[(y & 1) * 2 + (x & 1)]; + let src = img.data.clone(); + for &i in bad { + let (x, y) = (i % w, i / w); + let me = color(x, y); + for radius in [2isize, 4, 6, 8] { + let (mut sum, mut n) = (0f64, 0u32); + for dy in -radius..=radius { + for dx in -radius..=radius { + let (sx, sy) = (x as isize + dx, y as isize + dy); + if sx < 0 || sy < 0 || sx as usize >= w || sy as usize >= h || (dx == 0 && dy == 0) { + continue; + } + let (sx, sy) = (sx as usize, sy as usize); + if color(sx, sy) == me && !is_bad[sy * w + sx] { + sum += f64::from(src[sy * w + sx]); + n += 1; + } + } + } + if n > 0 { + if let Some(v) = img.data.get_mut(i) { + *v = (sum / f64::from(n)) as f32; + } + break; + } + } + } + Ok(()) +} + +fn trim(r: &mut Reader, img: &mut Image) -> Result<()> { + let (t, l, b, rt) = (r.u32()? as usize, r.u32()? as usize, r.u32()? as usize, r.u32()? as usize); + if !img.allow_trim { + return Err(Error::Unsupported("TrimBounds before the demosaic".into())); + } + let (b, rt) = (b.min(img.height), rt.min(img.width)); + if t >= b || l >= rt { + return Err(Error::Malformed("TrimBounds outside the picture")); + } + let (np, w) = (img.planes, img.width); + let nw = rt - l; + let mut out = Vec::with_capacity(nw * (b - t) * np); + for y in t..b { + out.extend_from_slice(img.data.get((y * w + l) * np..(y * w + rt) * np).ok_or(Error::Malformed("TrimBounds"))?); + } + *img.data = out; + img.width = nw; + img.height = b - t; + img.origin = (img.origin.0 + l, img.origin.1 + t); + Ok(()) +} + +/// Catmull–Rom weights for a fractional position `t` in 0–1. +fn cubic(t: f64) -> [f64; 4] { + let (t2, t3) = (t * t, t * t * t); + [(-t3 + 2.0 * t2 - t) / 2.0, (3.0 * t3 - 5.0 * t2 + 2.0) / 2.0, (-3.0 * t3 + 4.0 * t2 + t) / 2.0, (t3 - t2) / 2.0] +} + +/// The radial model of a warp opcode for one plane. +enum Radial { + /// WarpRectilinear: f = kr0 + kr1 r² + kr2 r⁴ + kr3 r⁶, tangential kt0, kt1. + Rectilinear([f64; 4], [f64; 2]), + /// WarpFisheye: f = (kr0 θ + kr1 θ³ + kr2 θ⁵ + kr3 θ⁷) / r, θ = atan r. + Fisheye([f64; 4]), + /// WarpRectilinear2: f(r) = Σ kr_p r^p (p ≤ 14) at r clamped to the valid range, multiplied + /// (or divided, reciprocal) into r; tangential kt0, kt1. + Rectilinear2 { kr: [f64; 15], kt: [f64; 2], min: f64, max: f64, reciprocal: bool }, +} + +impl Radial { + /// Source offset (normalised) of destination offset (dx, dy). + fn map(&self, dx: f64, dy: f64) -> (f64, f64) { + let r2 = dx * dx + dy * dy; + let tangential = |kt: &[f64; 2]| (kt[0] * (2.0 * dx * dy) + kt[1] * (r2 + 2.0 * dx * dx), kt[1] * (2.0 * dx * dy) + kt[0] * (r2 + 2.0 * dy * dy)); + match self { + Radial::Rectilinear(kr, kt) => { + let f = kr[0] + r2 * (kr[1] + r2 * (kr[2] + r2 * kr[3])); + let (tx, ty) = tangential(kt); + (f * dx + tx, f * dy + ty) + } + Radial::Fisheye(kr) => { + let r = r2.sqrt(); + let f = if r < 1e-12 { + kr[0] + } else { + let t = r.atan(); + let t2 = t * t; + t * (kr[0] + t2 * (kr[1] + t2 * (kr[2] + t2 * kr[3]))) / r + }; + (f * dx, f * dy) + } + Radial::Rectilinear2 { kr, kt, min, max, reciprocal } => { + let r = r2.sqrt().clamp(*min, *max); + let f = kr.iter().rev().fold(0.0, |acc, k| acc * r + k); + let f = if *reciprocal { + if f.abs() < 1e-12 { + return (dx, dy); + } + 1.0 / f + } else { + f + }; + let (tx, ty) = tangential(kt); + (f * dx + tx, f * dy + ty) + } + } + } +} + +fn warp(id: u32, r: &mut Reader, img: &mut Image) -> Result<()> { + if img.mosaic { + return Err(Error::Unsupported("a warp opcode on a CFA mosaic".into())); + } + let per_set = match id { + 1 => 6, + 2 => 4, + _ => 19, + }; + let n = r.count(per_set * 8)?; + if n != 1 && n != img.planes { + return Err(Error::Malformed("warp opcode coefficient sets")); + } + let mut sets = Vec::with_capacity(n); + for _ in 0..n { + let v: Vec = (0..per_set).map(|_| r.f64()).collect::>()?; + sets.push(match id { + 1 => Radial::Rectilinear([v[0], v[1], v[2], v[3]], [v[4], v[5]]), + 2 => Radial::Fisheye([v[0], v[1], v[2], v[3]]), + _ => { + let mut kr = [0f64; 15]; + kr.copy_from_slice(&v[..15]); + let (min, max) = (v[17].clamp(0.0, 1.0), v[18].clamp(0.0, 1.0)); + if min >= max { + return Err(Error::Malformed("WarpRectilinear2 valid radius range")); + } + Radial::Rectilinear2 { kr, kt: [v[15], v[16]], min, max, reciprocal: false } + } + }); + } + let (cx, cy) = (r.f64()?, r.f64()?); + if id == 14 { + let reciprocal = r.u32()? != 0; + for s in &mut sets { + if let Radial::Rectilinear2 { reciprocal: rr, .. } = s { + *rr = reciprocal; + } + } + } + let (cx, cy, m) = centre(img, cx, cy); + if m <= 0.0 { + return Ok(()); + } + let (w, h, np) = (img.width, img.height, img.planes); + let src = img.data.clone(); + let sample = |p: usize, x: isize, y: isize| { + let (x, y) = (x.clamp(0, w as isize - 1) as usize, y.clamp(0, h as isize - 1) as usize); + f64::from(src.get((y * w + x) * np + p).copied().unwrap_or(0.0)) + }; + for (i, px) in img.data.chunks_exact_mut(np.max(1)).enumerate() { + let (x, y) = ((i % w) as f64, (i / w) as f64); + let (dx, dy) = ((x - cx) / m, (y - cy) / m); + for (p, v) in px.iter_mut().enumerate() { + let model = sets.get(p).or(sets.first()); + let Some(model) = model else { continue }; + let (ox, oy) = model.map(dx, dy); + let (sx, sy) = (cx + m * ox, cy + m * oy); + if !sx.is_finite() || !sy.is_finite() { + continue; + } + let (fx, fy) = (sx.floor(), sy.floor()); + let (wx, wy) = (cubic(sx - fx), cubic(sy - fy)); + let (ix, iy) = (fx.clamp(-4.0, w as f64 + 4.0) as isize, fy.clamp(-4.0, h as f64 + 4.0) as isize); + let mut acc = 0.0; + for (j, wyj) in wy.iter().enumerate() { + let mut row = 0.0; + for (k, wxk) in wx.iter().enumerate() { + row += wxk * sample(p, ix + k as isize - 1, iy + j as isize - 1); + } + acc += wyj * row; + } + *v = acc as f32; + } + } + Ok(()) +} + +#[cfg(any(test, feature = "writer"))] +/// Encode an opcode list (big-endian), for tests: (id, flags, parameters). +pub fn encode(list: &[(u32, u32, Vec)]) -> Vec { + let mut out = (list.len() as u32).to_be_bytes().to_vec(); + for (id, flags, params) in list { + let version: u32 = if *id == 14 { 0x0106_0000 } else { 0x0103_0000 }; + for v in [*id, version, *flags, params.len() as u32] { + out.extend(v.to_be_bytes()); + } + out.extend(params); + } + out +} + +#[cfg(test)] +mod tests { + use super::*; + + fn u(v: &[u32]) -> Vec { + v.iter().flat_map(|x| x.to_be_bytes()).collect() + } + + fn d(v: &[f64]) -> Vec { + v.iter().flat_map(|x| x.to_be_bytes()).collect() + } + + fn f(v: &[f32]) -> Vec { + v.iter().flat_map(|x| x.to_be_bytes()).collect() + } + + fn area(rect_plane_pitch: [u32; 8]) -> Vec { + u(&rect_plane_pitch) + } + + /// Run one opcode on a float picture (range 0–1). + fn run1(id: u32, params: Vec, data: &mut Vec, w: usize, h: usize, planes: usize) -> Result<(usize, usize)> { + let list = parse(&encode(&[(id, 0, params)]))?; + let mut img = Image { data, width: w, height: h, planes, max: 1.0, mosaic: false, allow_trim: true, origin: (0, 0) }; + apply(&list, &mut img, false)?; + Ok((img.width, img.height)) + } + + fn close(a: f32, b: f32) -> bool { + (a - b).abs() < 1e-5 + } + + #[test] + fn list_parsing_and_hostile_lists() { + let l = encode(&[(9, 1, vec![1, 2, 3]), (3, 2, vec![])]); + let p = parse(&l).unwrap(); + assert_eq!(p.len(), 2); + assert_eq!((p[0].id, p[0].version, p[0].params.as_slice()), (9, 0x0103_0000, &[1u8, 2, 3][..])); + assert!(p[0].optional() && !p[0].preview_skippable()); + assert!(!p[1].optional() && p[1].preview_skippable()); + // every truncation is an error, never a panic + for n in 0..l.len() { + assert!(parse(&l[..n]).is_err(), "{n}"); + } + // counts and sizes beyond the data + assert!(parse(&u(&[u32::MAX])).is_err()); + assert!(parse(&u(&[1, 9, 0, 0, u32::MAX])).is_err()); + assert!(parse(&u(&[(MAX_OPCODES + 1) as u32]).into_iter().chain(vec![0; 16 * (MAX_OPCODES + 1)]).collect::>()).is_err()); + assert_eq!(parse(&u(&[0])).unwrap(), vec![]); + } + + #[test] + fn unknown_and_newer_opcodes() { + let mut data = vec![0.5f32; 4]; + let mut img = Image { data: &mut data, width: 2, height: 2, planes: 1, max: 1.0, mosaic: false, allow_trim: false, origin: (0, 0) }; + // unknown id: skipped when optional, an error when not + assert!(apply(&parse(&encode(&[(99, 1, vec![])])).unwrap(), &mut img, false).is_ok()); + assert!(matches!(apply(&parse(&encode(&[(99, 0, vec![])])).unwrap(), &mut img, false), Err(Error::Unsupported(_)))); + // a known id from a newer DNG version is unknown + let newer = [Opcode { id: 3, version: 0x0108_0000, flags: 0, params: vec![] }]; + assert!(matches!(apply(&newer, &mut img, false), Err(Error::Unsupported(_)))); + // an opcode this stage cannot run: skipped when optional + let trim = Opcode { id: 6, version: 0x0103_0000, flags: 1, params: u(&[0, 0, 1, 1]) }; + assert!(apply(std::slice::from_ref(&trim), &mut img, false).is_ok()); + assert!(apply(&[Opcode { flags: 0, ..trim }], &mut img, false).is_err()); + // preview quality skips flagged opcodes + let scale = Opcode { id: 12, version: 0x0103_0000, flags: 2, params: [area([0, 0, 2, 2, 0, 1, 1, 1]), u(&[2]), f(&[0.5, 0.5])].concat() }; + apply(std::slice::from_ref(&scale), &mut img, true).unwrap(); + assert_eq!(img.data[0], 0.5); + apply(&[scale], &mut img, false).unwrap(); + assert_eq!(img.data[0], 0.25); + // malformed parameters are errors + assert!(apply(&[Opcode { id: 10, version: 0x0103_0000, flags: 0, params: vec![0; 5] }], &mut img, false).is_err()); + } + + #[test] + fn map_table_polynomial_delta_and_scale() { + // 4 × 2, two planes; affect plane 1 only, every second column + let mut data: Vec = (0..16).map(|i| i as f32 / 20.0).collect(); + let orig = data.clone(); + // MapTable: 16-bit identity reversed (v → 1 − v) + let table: Vec = (0..=65_535u32).rev().flat_map(|v| (v as u16).to_be_bytes()).collect(); + run1(7, [area([0, 0, 2, 4, 1, 1, 1, 2]), u(&[65_536]), table].concat(), &mut data, 4, 2, 2).unwrap(); + for (i, (a, b)) in data.iter().zip(&orig).enumerate() { + let (x, p) = ((i / 2) % 4, i % 2); + let want = if p == 1 && x % 2 == 0 { 1.0 - (b * 65_535.0).round() / 65_535.0 } else { *b }; + assert!(close(*a, want), "{i}: {a} {want}"); + } + // MapPolynomial 0.1 + 0.5 x on everything + let mut data = orig.clone(); + run1(8, [area([0, 0, 2, 4, 0, 2, 1, 1]), u(&[1]), d(&[0.1, 0.5])].concat(), &mut data, 4, 2, 2).unwrap(); + for (a, b) in data.iter().zip(&orig) { + assert!(close(*a, 0.1 + 0.5 * b)); + } + // DeltaPerRow (rows 0, 1), ScalePerColumn (columns 1..3) + let mut data = orig.clone(); + run1(10, [area([0, 0, 2, 4, 0, 2, 1, 1]), u(&[2]), f(&[0.01, -0.02])].concat(), &mut data, 4, 2, 2).unwrap(); + for (i, (a, b)) in data.iter().zip(&orig).enumerate() { + let want = (b + if i < 8 { 0.01 } else { -0.02 }).clamp(0.0, 1.0); + assert!(close(*a, want), "{i}"); + } + let mut data = orig.clone(); + run1(13, [area([0, 1, 2, 3, 0, 2, 1, 1]), u(&[2]), f(&[2.0, 0.5])].concat(), &mut data, 4, 2, 2).unwrap(); + for (i, (a, b)) in data.iter().zip(&orig).enumerate() { + let x = (i / 2) % 4; + let want = match x { + 1 => (b * 2.0).min(1.0), + 2 => b * 0.5, + _ => *b, + }; + assert!(close(*a, want), "{i}"); + } + // too few deltas for the rows is an error + let mut data = orig.clone(); + assert!(run1(10, [area([0, 0, 2, 4, 0, 2, 1, 1]), u(&[1]), f(&[0.1])].concat(), &mut data, 4, 2, 2).is_err()); + // zero pitch is an error, areas outside the picture are clamped + assert!(run1(12, [area([0, 0, 2, 4, 0, 1, 0, 1]), u(&[0])].concat(), &mut data, 4, 2, 2).is_err()); + let mut data = orig.clone(); + run1(12, [area([1, 3, 99, 99, 5, 9, 1, 1]), u(&[1]), f(&[0.0])].concat(), &mut data, 4, 2, 2).unwrap(); + assert_eq!(data, orig, "plane 5 does not exist"); + } + + #[test] + fn gain_map_interpolates_and_replicates() { + // 2 × 2 map at the corners of the picture's relative coordinates 0 and 0.5 + let (w, h) = (8, 8); + let mut data = vec![0.25f32; w * h]; + let gains = f(&[1.0, 2.0, 3.0, 4.0]); + run1(9, [area([0, 0, 8, 8, 0, 1, 1, 1]), u(&[2, 2]), d(&[0.5, 0.5, 0.0, 0.0]), u(&[1]), gains].concat(), &mut data, w, h, 1).unwrap(); + let g = |x: usize, y: usize| data[y * w + x] / 0.25; + assert!(close(g(0, 0), 1.0)); + assert!(close(g(4, 0), 2.0), "{}", g(4, 0)); + assert!(close(g(0, 4), 3.0)); + assert!(close(g(2, 2), 2.5), "bilinear in the middle: {}", g(2, 2)); + assert!(close(g(7, 7), 4.0), "replicated past the map"); + // a map larger than its parameters is an error + let mut data = vec![0.25f32; 4]; + assert!( + run1(9, [area([0, 0, 2, 2, 0, 1, 1, 1]), u(&[1000, 1000]), d(&[0.1, 0.1, 0.0, 0.0]), u(&[1]), f(&[1.0])].concat(), &mut data, 2, 2, 1).is_err() + ); + } + + #[test] + fn vignette_and_warps() { + // FixVignetteRadial k0 = 1: g = 1 + r², the centre unchanged, the corners doubled + let (w, h) = (5, 5); + let mut data = vec![0.25f32; w * h]; + run1(3, d(&[1.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5]), &mut data, w, h, 1).unwrap(); + assert!(close(data[2 * w + 2], 0.25)); + assert!(close(data[0], 0.5) && close(data[w * h - 1], 0.5)); + assert!(close(data[2 * w], 0.25 * 1.5), "r² = 0.5 at the edge middle: {}", data[2 * w]); + // identity warps leave the picture alone; kr0 = 0.5 magnifies a ramp about the centre + let (w, h) = (9, 7); + let ramp: Vec = (0..w * h).flat_map(|i| [(i % w) as f32 / 10.0, (i / w) as f32 / 10.0]).collect(); + let mut data = ramp.clone(); + run1(1, [u(&[1]), d(&[1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5])].concat(), &mut data, w, h, 2).unwrap(); + assert!(data.iter().zip(&ramp).all(|(a, b)| close(*a, *b))); + let mut data = ramp.clone(); + run1(1, [u(&[1]), d(&[0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5])].concat(), &mut data, w, h, 2).unwrap(); + // x' = 4 + 0.5 (x − 4): the cubic reproduces a linear ramp exactly + for x in 1..w - 1 { + let want = (4.0 + 0.5 * (x as f32 - 4.0)) / 10.0; + assert!(close(data[(3 * w + x) * 2], want), "x {x}: {} vs {want}", data[(3 * w + x) * 2]); + } + // one coefficient set per plane: plane 1 untouched + let mut data = ramp.clone(); + run1(1, [u(&[2]), d(&[0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5])].concat(), &mut data, w, h, 2).unwrap(); + assert!(data.iter().skip(1).step_by(2).zip(ramp.iter().skip(1).step_by(2)).all(|(a, b)| close(*a, *b))); + // three sets for two planes is malformed + let mut data = ramp.clone(); + assert!(run1(1, [u(&[3]), vec![0; 3 * 48 + 16]].concat(), &mut data, w, h, 2).is_err()); + // WarpFisheye kr0 = 1: f = atan(r) / r pulls sources towards the centre + let mut data = ramp.clone(); + run1(2, [u(&[1]), d(&[1.0, 0.0, 0.0, 0.0, 0.5, 0.5])].concat(), &mut data, w, h, 2).unwrap(); + let m = (16.0f64 + 9.0).sqrt(); + let dx = 3.0 / m; + let want = (4.0 + m * dx.atan()) / 10.0; + assert!((f64::from(data[(3 * w + 7) * 2]) - want).abs() < 1e-3, "{} vs {want}", data[(3 * w + 7) * 2]); + // WarpRectilinear2 identity (kr0 = 1), and the skip rule for a following warp + let wr2 = |flags| (14, flags, [u(&[1]), d(&[1.0]), d(&[0.0; 14]), d(&[0.0, 0.0, 0.0, 1.0, 0.5, 0.5]), u(&[0])].concat()); + let halve = (1, 0, [u(&[1]), d(&[0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5])].concat()); + for (flags, skipped) in [(1, true), (0, false)] { + let list = parse(&encode(&[wr2(flags), halve.clone()])).unwrap(); + let mut data = ramp.clone(); + let mut img = Image { data: &mut data, width: w, height: h, planes: 2, max: 1.0, mosaic: false, allow_trim: false, origin: (0, 0) }; + apply(&list, &mut img, false).unwrap(); + let unchanged = data.iter().zip(&ramp).all(|(a, b)| close(*a, *b)); + assert_eq!(unchanged, skipped, "optional WarpRectilinear2 {flags}"); + } + // warps are refused on a mosaic + let mut data = vec![0.5f32; 16]; + let mut img = Image { data: &mut data, width: 4, height: 4, planes: 1, max: 1.0, mosaic: true, allow_trim: false, origin: (0, 0) }; + assert!(apply(&parse(&encode(&[halve])).unwrap(), &mut img, false).is_err()); + } + + #[test] + fn bad_pixels_and_trim() { + // RGGB 6 × 6 with each colour flat; two bad pixels and a bad 2 × 2 rectangle + let (w, h) = (6, 6); + let level = [0.2f32, 0.5, 0.5, 0.8]; + let flat: Vec = (0..w * h).map(|i| level[(i / w % 2) * 2 + i % 2]).collect(); + let mut data = flat.clone(); + data[2 * w + 2] = 0.0; + data[3 * w + 4] = 0.0; + // FixBadPixelsConstant 0, phase 0 (red top left), integer stage + let mut ints: Vec = data.iter().map(|v| (v * 65_535.0).round()).collect(); + let list = parse(&encode(&[(4, 0, u(&[0, 0]))])).unwrap(); + let mut img = Image { data: &mut ints, width: w, height: h, planes: 1, max: 65_535.0, mosaic: true, allow_trim: false, origin: (0, 0) }; + apply(&list, &mut img, false).unwrap(); + for (a, b) in ints.iter().zip(&flat) { + assert_eq!(*a, (b * 65_535.0).round()); + } + // FixBadPixelsList: a point and a rectangle + let mut data = flat.clone(); + for i in [w + 1, 4 * w + 3, 4 * w + 4, 5 * w + 3, 5 * w + 4] { + data[i] = 0.9; + } + run1(5, u(&[0, 1, 1, 1, 1, 4, 3, 6, 5]), &mut data, w, h, 1).unwrap(); + assert!(data.iter().zip(&flat).all(|(a, b)| close(*a, *b)), "{data:?}"); + // a bad phase is malformed; huge lists are refused before allocating + let mut data = flat.clone(); + assert!(run1(5, u(&[7, 1, 0, 1, 1]), &mut data, w, h, 1).is_err()); + assert!(run1(5, u(&[0, u32::MAX, u32::MAX]), &mut data, w, h, 1).is_err()); + // TrimBounds + let mut data: Vec = (0..w * h * 3).map(|i| i as f32 / 200.0).collect(); + let (tw, th) = run1(6, u(&[1, 2, 4, 5]), &mut data, w, h, 3).unwrap(); + assert_eq!((tw, th, data.len()), (3, 3, 27)); + assert!(close(data[0], ((w + 2) * 3) as f32 / 200.0)); + let mut data = vec![0f32; 4]; + assert!(run1(6, u(&[1, 1, 1, 2]), &mut data, 2, 2, 1).is_err()); + } +} diff --git a/crates/dng/src/options.rs b/crates/dng/src/options.rs new file mode 100644 index 0000000000..20004c3ff8 --- /dev/null +++ b/crates/dng/src/options.rs @@ -0,0 +1,84 @@ +//! Raw development settings ([`ProcessOptions`]): what a "Raw settings" panel can change about how +//! a DNG is developed. `ProcessOptions::default()` develops exactly as the camera file says (the +//! as-shot white balance, no exposure change, highlights clipped, the high-quality demosaic, +//! scene-linear Rec. 709, full resolution). Feature `serde` derives Serialize / Deserialize (missing +//! fields take their defaults). + +/// White balance. +#[derive(Clone, Copy, Debug, Default, PartialEq)] +#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))] +pub enum WhiteBalance { + /// The camera's AsShotNeutral / AsShotWhiteXY (D65 without either). + #[default] + AsShot, + /// A white of `temperature` kelvin (1 000–15 000, on the Planckian locus) and `tint` + /// (positive makes the picture more magenta, negative greener, as in common raw converters; + /// ±150 is the useful range, D65 is about +10, see [`crate::color::TINT_PER_DUV`]). Needs the file's ColorMatrix; a DNG without one has no + /// colour calibration to place the white in and keeps its as-shot balance. + Custom { temperature: f64, tint: f64 }, +} + +/// What happens to values the sensor clipped. +#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)] +#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))] +pub enum Highlights { + /// White-balanced samples are clipped at the sensor's white, so clipped highlights stay + /// neutral (no magenta or cyan casts). + #[default] + Clip, + /// Reconstruct clipped channels from the unclipped ones. **Not implemented yet: develops + /// exactly like [`Highlights::Clip`]**, so the setting can already be stored and offered. + Recover, +} + +/// Demosaic quality (CFA images only). +#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)] +#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))] +pub enum DemosaicQuality { + /// Bilinear interpolation: softer, a little faster (playback, scrubbing). + Fast, + /// Malvar–He–Cutler gradient-corrected interpolation for 2 × 2 Bayer patterns. + #[default] + High, +} + +/// The colour space of developed pictures. +#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)] +#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))] +pub enum OutputSpace { + /// Scene-linear light with Rec. 709 / sRGB primaries and a D65 white; a white surface the + /// sensor saw at its clipping level is 1.0 × 2^(BaselineExposure + exposure). + #[default] + SceneLinearRec709, +} + +/// How to develop a DNG ([`crate::Dng::develop_with`]). +#[derive(Clone, Copy, Debug, Default, PartialEq)] +#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))] +#[cfg_attr(feature = "serde", serde(default))] +pub struct ProcessOptions { + pub white_balance: WhiteBalance, + /// Exposure change in stops, added to the file's BaselineExposure (the sum is limited to + /// ±16 stops; a non-finite value counts as 0). + pub exposure: f64, + pub highlights: Highlights, + pub demosaic: DemosaicQuality, + pub output: OutputSpace, + /// The smallest picture (width, height) the caller needs, e.g. a monitor or a proxy size. + /// `None` develops at full resolution. With a size, the decoder may develop a smaller + /// picture (half resolution when both dimensions still cover the target) that is never + /// smaller than the target; [`crate::Image`]'s `width` and `height` give the size made, + /// [`crate::Info`] still the full one. + pub target_size: Option<(usize, usize)>, + /// Apply the file's OpcodeList1–3 (lens shading, warps, bad pixels, trims). **Off by default**: + /// the opcode interpretations have not yet been compared with real camera files (see the + /// README), so a file's lists are ignored unless this is set. + pub apply_opcodes: bool, +} + +impl ProcessOptions { + /// The exposure change in stops (0 for a non-finite value). + pub fn exposure_stops(&self) -> f64 { + if self.exposure.is_finite() { self.exposure } else { 0.0 } + } +} diff --git a/crates/dng/src/process.rs b/crates/dng/src/process.rs index 41d9b25aed..fc288babf5 100644 --- a/crates/dng/src/process.rs +++ b/crates/dng/src/process.rs @@ -18,7 +18,9 @@ use rayon::prelude::*; use crate::color::{M3, V3}; -use crate::demosaic::{Pattern, demosaic}; +use crate::demosaic::Pattern; +use crate::opcode::{self, Opcode}; +use crate::options::ProcessOptions; use crate::raw::RawImage; use crate::{Error, Result}; @@ -40,6 +42,8 @@ pub struct Params { pub crop_size: Option<(f64, f64)>, /// CFA pattern in camera planes (`None`: LinearRaw). pub cfa: Option, + /// OpcodeList1–3. + pub opcodes: crate::opcode::Lists, } /// A rectangle: x, y, width, height. @@ -70,8 +74,69 @@ pub fn crop_rect(p: &Params, aw: usize, ah: usize) -> Rect { (x, y, w, h) } -/// Normalised, white-balanced samples of the active area (steps 1–4), one plane per sample. -fn normalise(raw: &RawImage, p: &Params, (ax, ay, aw, ah): Rect, neutral: &V3, planes: usize) -> Result> { +/// Rows-parallel map over `out` in rows of `row_len` values (sequential without `threads`). +fn for_rows(out: &mut [f32], row_len: usize, f: impl Fn(usize, &mut [f32]) + Sync + Send) { + let row_len = row_len.max(1); + #[cfg(feature = "threads")] + out.par_chunks_mut(row_len).enumerate().for_each(|(y, o)| f(y, o)); + #[cfg(not(feature = "threads"))] + out.chunks_mut(row_len).enumerate().for_each(|(y, o)| f(y, o)); +} + +fn gcd(a: usize, b: usize) -> usize { + if b == 0 { a } else { gcd(b, a % b) } +} + +fn lcm(a: usize, b: usize) -> usize { + a / gcd(a, b).max(1) * b +} + +/// Most values in a normalisation lookup table (sites × codes). +const MAX_LUT: usize = 1 << 21; + +/// Steps 2–4 for one sample: linearisation, black and white levels, and (with `wb`) white +/// balance with the clip at 1. +struct Normaliser<'a> { + spp: usize, + br: usize, + bc: usize, + black: Vec, + white: Vec, + lin: Option<&'a [u16]>, + inv_neutral: Vec, + cfa: Option<&'a Pattern>, + wb: bool, +} + +impl Normaliser<'_> { + #[inline] + fn value(&self, code: u16, x: usize, y: usize, s: usize, dh: f64, dv: f64) -> f32 { + let v = match self.lin { + Some(t) => f64::from(t.get(usize::from(code)).or(t.last()).copied().unwrap_or(0)), + None => f64::from(code), + }; + let b = self.black.get(((y % self.br) * self.bc + x % self.bc) * self.spp + s).copied().unwrap_or(0.0) + dh + dv; + let range = self.white.get(s).copied().unwrap_or(1.0) - b; + let n = if range > 0.0 { ((v - b) / range).clamp(0.0, 1.0) } else { 0.0 }; + if !self.wb { + return n as f32; + } + let plane = match self.cfa { + Some(c) => usize::from(c.color(x, y)), + None => s, + }; + let wb = self.inv_neutral.get(plane).copied().unwrap_or(1.0); + (n as f32 * wb).min(1.0) + } +} + +/// Normalised samples of the active area (steps 1–4), one plane per sample; white balanced and +/// clipped at 1 when `wb` (else only clipped to 0–1, for OpcodeList2). +/// +/// Without black-level deltas every sample is a function of its code and its position in the +/// black-level and CFA repeat, so a lookup table per site replaces the arithmetic (same values: +/// each entry is computed by the same expression). +fn normalise_with(raw: &RawImage, p: &Params, (ax, ay, aw, ah): Rect, neutral: &V3, planes: usize, wb: bool) -> Result> { let spp = raw.samples_per_pixel; let (br, bc) = (p.black_repeat.0.max(1), p.black_repeat.1.max(1)); if br > 64 || bc > 64 { @@ -84,44 +149,73 @@ fn normalise(raw: &RawImage, p: &Params, (ax, ay, aw, ah): Rect, neutral: &V3, p n if n >= need => p.black[..need].to_vec(), _ => return Err(Error::Malformed("BlackLevel has too few values")), }; - let max_code = ((1u32 << raw.bits) - 1) as f64; + let max_code = ((1u32 << raw.bits.min(16)) - 1) as f64; let white: Vec = (0..spp).map(|s| p.white.get(s).or(p.white.first()).copied().unwrap_or(max_code)).collect(); let lin = p.linearization.as_deref().filter(|t| !t.is_empty()); let inv_neutral: Vec = (0..planes).map(|c| (1.0 / neutral.get(c).copied().unwrap_or(1.0)) as f32).collect(); let cfa = p.cfa.as_ref(); + let n = Normaliser { spp, br, bc, black, white, lin, inv_neutral, cfa, wb }; + let row_len = aw * spp; + let src = |y: usize| raw.data.get((ay + y) * raw.width * spp + ax * spp..).unwrap_or(&[]); let mut out = vec![0f32; aw * ah * spp]; - let row = |y: usize, o: &mut [f32]| { - let dv = p.black_delta_v.get(y).copied().unwrap_or(0.0); - let src_row = (ay + y) * raw.width * spp; - for x in 0..aw { - let dh = p.black_delta_h.get(x).copied().unwrap_or(0.0); - for s in 0..spp { - let code = raw.data.get(src_row + (ax + x) * spp + s).copied().unwrap_or(0); - let v = match lin { - Some(t) => f64::from(t.get(usize::from(code)).or(t.last()).copied().unwrap_or(0)), - None => f64::from(code), - }; - let b = black[((y % br) * bc + x % bc) * spp + s] + dh + dv; - let range = white[s] - b; - let n = if range > 0.0 { ((v - b) / range).clamp(0.0, 1.0) } else { 0.0 }; - let plane = match cfa { - Some(c) => usize::from(c.color(x, y)), - None => s, + // the repeat of the black levels and the CFA together + let (rr, rc) = match cfa { + Some(c) => (lcm(br, c.rows.max(1)), lcm(bc, c.cols.max(1))), + None => (br, bc), + }; + let deltas = p.black_delta_h.iter().chain(&p.black_delta_v).any(|&d| d != 0.0); + let lut_len = 1usize << raw.bits.clamp(1, 16); + let sites = rr * rc * spp; + let use_lut = !deltas && sites.saturating_mul(lut_len) <= MAX_LUT && sites * lut_len <= out.len(); + if use_lut { + // lut[((y % rr) * rc + x % rc) * spp + s][code] + let mut lut = vec![0f32; sites * lut_len]; + let fill = |site: usize, t: &mut [f32]| { + let (s, xy) = (site % spp, site / spp); + let (y, x) = (xy / rc, xy % rc); + for (code, v) in t.iter_mut().enumerate() { + *v = n.value(code as u16, x, y, s, 0.0, 0.0); + } + }; + for_rows(&mut lut, lut_len, fill); + let lut = &lut; + for_rows(&mut out, row_len, |y, o| { + let src = src(y); + let row_sites = (y % rr) * rc * spp; + let mut site = row_sites; + let row_end = row_sites + rc * spp; + for (x, (o, &code)) in o.iter_mut().zip(src).enumerate() { + let table = lut.get(site * lut_len..(site + 1) * lut_len).unwrap_or(&[]); + *o = match table.get(usize::from(code)) { + Some(&v) => v, + None => n.value(code, x / spp, y, x % spp, 0.0, 0.0), }; - let wb = inv_neutral.get(plane).copied().unwrap_or(1.0); - if let Some(o) = o.get_mut(x * spp + s) { - *o = (n as f32 * wb).min(1.0); + site += 1; + if site == row_end { + site = row_sites; } } - } - }; - #[cfg(feature = "threads")] - out.par_chunks_mut((aw * spp).max(1)).enumerate().for_each(|(y, o)| row(y, o)); - #[cfg(not(feature = "threads"))] - out.chunks_mut((aw * spp).max(1)).enumerate().for_each(|(y, o)| row(y, o)); + }); + } else { + for_rows(&mut out, row_len, |y, o| { + let src = src(y); + let dv = p.black_delta_v.get(y).copied().unwrap_or(0.0); + for (x, (px, codes)) in o.chunks_exact_mut(spp).zip(src.chunks_exact(spp)).enumerate() { + let dh = p.black_delta_h.get(x).copied().unwrap_or(0.0); + for (s, (o, &code)) in px.iter_mut().zip(codes).enumerate() { + *o = n.value(code, x, y, s, dh, dv); + } + } + }); + } Ok(out) } +/// Normalised, white-balanced samples of the active area (steps 1–4), one plane per sample. +fn normalise(raw: &RawImage, p: &Params, rect: Rect, neutral: &V3, planes: usize) -> Result> { + normalise_with(raw, p, rect, neutral, planes, true) +} + /// A developed picture: premultiplied linear Rec. 709 RGBA floats (alpha 1). pub struct Developed { pub width: usize, @@ -131,6 +225,71 @@ pub struct Developed { /// Develop `raw` (steps 1–7). `neutral` and `matrix` come from [`crate::color::Calibration::solve`]. pub fn develop(raw: &RawImage, p: &Params, neutral: &V3, matrix: &M3) -> Result { + develop_with(raw, p, neutral, matrix, &ProcessOptions::default()) +} + +/// The output pixel of white-balanced camera values `c` (step 7). +#[inline(always)] +fn to_output(m: &[[f32; 3]; 3], mono: bool, c: [f32; 3], d: &mut [f32; 4]) { + let px = if mono { + // grey: the matrix's luminance gain only + let g = c[0] * (m[1][0] + m[1][1] + m[1][2]); + [g, g, g] + } else { + [0, 1, 2].map(|i| m[i][0] * c[0] + m[i][1] * c[1] + m[i][2] * c[2]) + }; + *d = [px[0].max(0.0), px[1].max(0.0), px[2].max(0.0), 1.0]; +} + +/// The reduction factor (1, 2, 4 or 8) for a `cw` × `ch` picture when the caller needs at least +/// `target`: the largest that still covers it. +pub fn reduction(target: Option<(usize, usize)>, cw: usize, ch: usize) -> usize { + let Some((tw, th)) = target else { return 1 }; + let mut f = 1; + while f < 8 && cw / (f * 2) >= tw.max(1) && ch / (f * 2) >= th.max(1) { + f *= 2; + } + f +} + +/// Whether `list` has an opcode to apply (`preview`: those marked as needed only for full +/// quality are left out). +fn active(list: &[Opcode], preview: bool) -> bool { + list.iter().any(|o| !(preview && o.preview_skippable())) +} + +/// Whether this development is a reduced-resolution preview: a target size smaller than the +/// crop, and no OpcodeList3 work that has to be done at full size. +pub fn preview(p: &Params, o: &ProcessOptions, info: &crate::Info) -> bool { + let reduced = reduction(o.target_size, info.width, info.height) > 1 && p.cfa.as_ref().is_none_or(Pattern::is_bayer); + reduced && !active(&p.opcodes.list3, true) +} + +/// OpcodeList1 on the stored samples (integer range 0–65535). +pub fn apply_list1(raw: &mut RawImage, p: &Params, preview: bool) -> Result<()> { + if !active(&p.opcodes.list1, preview) { + return Ok(()); + } + let mut data: Vec = raw.data.iter().map(|&v| f32::from(v)).collect(); + let mut img = opcode::Image { + data: &mut data, + width: raw.width, + height: raw.height, + planes: raw.samples_per_pixel, + max: 65_535.0, + mosaic: p.cfa.is_some(), + allow_trim: false, + origin: (0, 0), + }; + opcode::apply(&p.opcodes.list1, &mut img, preview)?; + raw.data = data.iter().map(|&v| v.clamp(0.0, 65_535.0) as u16).collect(); + Ok(()) +} + +/// [`develop`] with raw settings. White balance and exposure are already in `neutral` and +/// `matrix` ([`crate::Dng::develop_with`] solves them); this applies the demosaic quality and the +/// reduced-resolution path. +pub fn develop_with(raw: &RawImage, p: &Params, neutral: &V3, matrix: &M3, o: &ProcessOptions) -> Result { let spp = raw.samples_per_pixel; if p.cfa.is_some() && spp != 1 { return Err(Error::Malformed("a CFA image with more than one sample per pixel")); @@ -138,41 +297,233 @@ pub fn develop(raw: &RawImage, p: &Params, neutral: &V3, matrix: &M3) -> Result< if p.cfa.is_none() && !matches!(spp, 1 | 3) { return Err(Error::Unsupported(format!("LinearRaw with {spp} samples per pixel"))); } - let active = active_rect(p, raw.width, raw.height)?; - let (_, _, aw, ah) = active; + if raw.data.len() < raw.width.saturating_mul(raw.height).saturating_mul(spp) { + return Err(Error::Truncated); + } + let active_r = active_rect(p, raw.width, raw.height)?; + let (_, _, aw, ah) = active_r; let planes = if p.cfa.is_some() { 3 } else { spp }; - let norm = normalise(raw, p, active, neutral, planes)?; - // three samples per pixel - let rgb = match (&p.cfa, spp) { - (Some(pat), _) => demosaic(&norm, aw, ah, pat), - (None, 3) => norm, - _ => norm.iter().flat_map(|&v| [v, v, v]).collect(), - }; - let (cx, cy, cw, ch) = crop_rect(p, aw, ah); + let crop = crop_rect(p, aw, ah); let m: [[f32; 3]; 3] = matrix.map(|r| r.map(|v| v as f32)); - let mono = p.cfa.is_none() && spp == 1; + let f = reduction(o.target_size, crop.2, crop.3); + let preview = f > 1 && p.cfa.as_ref().is_none_or(Pattern::is_bayer) && !active(&p.opcodes.list3, true); + let norm = if active(&p.opcodes.list2, preview) { + // OpcodeList2 works on the linear reference values, before the white balance + let mut norm = normalise_with(raw, p, active_r, neutral, planes, false)?; + let mut img = + opcode::Image { data: &mut norm, width: aw, height: ah, planes: spp, max: 1.0, mosaic: p.cfa.is_some(), allow_trim: false, origin: (0, 0) }; + opcode::apply(&p.opcodes.list2, &mut img, preview)?; + let inv: Vec = (0..planes).map(|c| (1.0 / neutral.get(c).copied().unwrap_or(1.0)) as f32).collect(); + for_rows(&mut norm, aw * spp, |y, row| { + for (i, v) in row.iter_mut().enumerate() { + let plane = match &p.cfa { + Some(c) => usize::from(c.color(i, y)), + None => i % spp, + }; + *v = (*v * inv.get(plane).copied().unwrap_or(1.0)).min(1.0); + } + }); + norm + } else { + normalise(raw, p, active_r, neutral, planes)? + }; + if preview { + return Ok(binned(&norm, aw, spp, p.cfa.as_ref(), crop, f, &m)); + } + if active(&p.opcodes.list3, false) { + return with_list3(norm, aw, ah, spp, p, crop, o, &m, neutral); + } + Ok(full(norm, (aw, ah, spp), p.cfa.as_ref(), crop, o, &m)) +} + +/// Steps 5–7 at full resolution: demosaic, crop and colour, row by row without intermediate +/// pictures for 2 × 2 Bayer patterns and LinearRaw. +fn full( + norm: Vec, + (aw, ah, spp): (usize, usize, usize), + cfa: Option<&Pattern>, + (cx, cy, cw, ch): Rect, + o: &ProcessOptions, + m: &[[f32; 3]; 3], +) -> Developed { + let mono = cfa.is_none() && spp == 1; let mut rgba = vec![0f32; cw * ch * 4]; - let row = |y: usize, o: &mut [f32]| { - let src = (cy + y) * aw + cx; - for x in 0..cw { - let Some(c) = rgb.get((src + x) * 3..(src + x) * 3 + 3) else { continue }; - let px = if mono { - // grey: the matrix's luminance gain only - let g = c[0] * (m[1][0] + m[1][1] + m[1][2]); - [g, g, g] + match cfa { + Some(pat) if pat.is_bayer() && aw >= 3 && ah >= 3 => { + let norm = &norm; + for_rows(&mut rgba, cw * 4, |y, out| { + let mut out = out.as_chunks_mut::<4>().0.iter_mut(); + crate::demosaic::BayerRow::new(norm, aw, ah, pat, o.demosaic, cy + y).span(cx, cx + cw, |rgb| { + if let Some(d) = out.next() { + to_output(m, false, rgb, d); + } + }); + }); + } + Some(pat) => { + let rgb = crate::demosaic::demosaic_with(&norm, aw, ah, pat, o.demosaic); + for_rows(&mut rgba, cw * 4, |y, out| { + let src = (cy + y) * aw + cx; + for (x, d) in out.as_chunks_mut::<4>().0.iter_mut().enumerate() { + if let Some(c) = rgb.get((src + x) * 3..(src + x) * 3 + 3) { + to_output(m, false, [c[0], c[1], c[2]], d); + } + } + }); + } + None => { + let norm = &norm; + for_rows(&mut rgba, cw * 4, |y, out| { + let src = (cy + y) * aw + cx; + for (x, d) in out.as_chunks_mut::<4>().0.iter_mut().enumerate() { + let c = if spp == 3 { + let Some(c) = norm.get((src + x) * 3..(src + x) * 3 + 3) else { continue }; + [c[0], c[1], c[2]] + } else { + let Some(&v) = norm.get(src + x) else { continue }; + [v, v, v] + }; + to_output(m, mono, c, d); + } + }); + } + } + Developed { width: cw, height: ch, rgba } +} + +/// Steps 5–7 with OpcodeList3: demosaic the whole active area, apply the list to the camera +/// values (white balance taken out again, range 0–1), then crop (within what TrimBounds left) +/// and colour. +#[allow(clippy::too_many_arguments)] +fn with_list3( + norm: Vec, + aw: usize, + ah: usize, + spp: usize, + p: &Params, + crop: Rect, + o: &ProcessOptions, + m: &[[f32; 3]; 3], + neutral: &V3, +) -> Result { + let mono = p.cfa.is_none() && spp == 1; + let planes = if mono { 1 } else { 3 }; + let mut rgb = match &p.cfa { + Some(pat) if pat.is_bayer() && aw >= 3 && ah >= 3 => { + let mut rgb = vec![0f32; aw * ah * 3]; + let norm = &norm; + for_rows(&mut rgb, aw * 3, |y, out| { + let mut out = out.as_chunks_mut::<3>().0.iter_mut(); + crate::demosaic::BayerRow::new(norm, aw, ah, pat, o.demosaic, y).span(0, aw, |c| { + if let Some(d) = out.next() { + *d = c; + } + }); + }); + rgb + } + Some(pat) => crate::demosaic::demosaic_with(&norm, aw, ah, pat, o.demosaic), + None => norm, + }; + let n: Vec = (0..planes).map(|c| if mono { 1.0 } else { neutral.get(c).copied().unwrap_or(1.0) as f32 }).collect(); + for px in rgb.chunks_exact_mut(planes) { + for (v, k) in px.iter_mut().zip(&n) { + *v *= k; + } + } + let mut img = opcode::Image { data: &mut rgb, width: aw, height: ah, planes, max: 1.0, mosaic: false, allow_trim: true, origin: (0, 0) }; + opcode::apply(&p.opcodes.list3, &mut img, false)?; + let (tw, th, (ox, oy)) = (img.width, img.height, img.origin); + for px in rgb.chunks_exact_mut(planes) { + for (v, k) in px.iter_mut().zip(&n) { + *v /= k; + } + } + // the default crop, inside the trimmed picture + let (cx, cy, cw, ch) = crop; + let (x0, y0) = (cx.max(ox), cy.max(oy)); + let (x1, y1) = ((cx + cw).min(ox + tw), (cy + ch).min(oy + th)); + if x0 >= x1 || y0 >= y1 { + return Err(Error::Malformed("TrimBounds outside the default crop")); + } + let (w, h) = (x1 - x0, y1 - y0); + let mut rgba = vec![0f32; w * h * 4]; + let rgb = &rgb; + for_rows(&mut rgba, w * 4, |y, out| { + let src = (y0 - oy + y) * tw + (x0 - ox); + for (x, d) in out.as_chunks_mut::<4>().0.iter_mut().enumerate() { + let i = (src + x) * planes; + let c = if mono { + let v = rgb.get(i).copied().unwrap_or(0.0); + [v, v, v] } else { - [0, 1, 2].map(|i| m[i][0] * c[0] + m[i][1] * c[1] + m[i][2] * c[2]) + let Some(c) = rgb.get(i..i + 3) else { continue }; + [c[0], c[1], c[2]] }; - if let Some(d) = o.get_mut(x * 4..x * 4 + 4) { - d.copy_from_slice(&[px[0].max(0.0), px[1].max(0.0), px[2].max(0.0), 1.0]); + to_output(m, mono, c, d); + } + }); + Ok(Developed { width: w, height: h, rgba }) +} + +/// Steps 5–7 at 1/`f` resolution: each output pixel averages every colour over an `f` × `f` +/// block of the crop (a 2 × 2 Bayer block has all three; no demosaic needed). +fn binned(norm: &[f32], aw: usize, spp: usize, cfa: Option<&Pattern>, (cx, cy, cw, ch): Rect, f: usize, m: &[[f32; 3]; 3]) -> Developed { + let mono = cfa.is_none() && spp == 1; + let (w, h) = ((cw / f).max(1), (ch / f).max(1)); + let mut rgba = vec![0f32; w * h * 4]; + // the colours of a 2 × 2 Bayer pattern (the only patterns binned), and how many samples of + // each an f × f block holds (f is even, so every block starts on the same phase) + let quad: [u8; 4] = cfa.map(|p| [p.color(0, 0), p.color(1, 0), p.color(0, 1), p.color(1, 1)]).unwrap_or([0; 4]); + let mut n = [0u32; 3]; + for by in 0..f { + for bx in 0..f { + match cfa { + Some(_) => n[usize::from(quad[((cy + by) & 1) * 2 + ((cx + bx) & 1)]).min(2)] += 1, + None => n = n.map(|v| v + 1), } } - }; - #[cfg(feature = "threads")] - rgba.par_chunks_mut(cw * 4).enumerate().for_each(|(y, o)| row(y, o)); - #[cfg(not(feature = "threads"))] - rgba.chunks_mut(cw * 4).enumerate().for_each(|(y, o)| row(y, o)); - Ok(Developed { width: cw, height: ch, rgba }) + } + let inv = n.map(|v| if v > 0 { 1.0 / v as f32 } else { 0.0 }); + for_rows(&mut rgba, w * 4, |y, out| { + let mut sums = vec![[0f32; 3]; w]; + for by in 0..f { + let ay = cy + y * f + by; + let start = (ay * aw + cx) * spp; + let Some(row) = norm.get(start..start + w * f * spp) else { continue }; + match cfa { + Some(_) => { + let even = usize::from(quad[(ay & 1) * 2 + (cx & 1)]).min(2); + let odd = usize::from(quad[(ay & 1) * 2 + ((cx + 1) & 1)]).min(2); + for (s, block) in sums.iter_mut().zip(row.chunks_exact(f)) { + for pair in block.as_chunks::<2>().0 { + s[even] += pair[0]; + s[odd] += pair[1]; + } + } + } + None if spp == 3 => { + for (s, block) in sums.iter_mut().zip(row.chunks_exact(f * 3)) { + for px in block.as_chunks::<3>().0 { + s[0] += px[0]; + s[1] += px[1]; + s[2] += px[2]; + } + } + } + None => { + for (s, block) in sums.iter_mut().zip(row.chunks_exact(f)) { + s[0] += block.iter().sum::(); + } + } + } + } + for (s, d) in sums.iter().zip(out.as_chunks_mut::<4>().0.iter_mut()) { + let c = if mono { [s[0] * inv[0]; 3] } else { [s[0] * inv[0], s[1] * inv[1], s[2] * inv[2]] }; + to_output(m, mono, c, d); + } + }); + Developed { width: w, height: h, rgba } } #[cfg(test)] diff --git a/crates/dng/src/raw.rs b/crates/dng/src/raw.rs index b69ba967e9..1ad135e19e 100644 --- a/crates/dng/src/raw.rs +++ b/crates/dng/src/raw.rs @@ -1,6 +1,9 @@ //! The raw image data of a DNG: finding the main raw IFD, and reading its strips or tiles //! (uncompressed or lossless JPEG) into one array of 16-bit samples. +#[cfg(feature = "threads")] +use rayon::prelude::*; + use crate::tiff::{Ifd, Tiff, tags}; use crate::{Error, Result, ljpeg}; @@ -21,7 +24,7 @@ pub struct RawLayout { pub height: usize, pub samples_per_pixel: usize, pub bits: u32, - /// 1 uncompressed, 7 lossless JPEG. + /// 1 uncompressed, 7 lossless JPEG, 34892 lossy JPEG. pub compression: u64, pub photometric: u64, } @@ -107,7 +110,12 @@ pub fn layout(t: &Tiff, ifd: &Ifd) -> Result { } 7 => {} 8 => return Err(Error::Unsupported("deflate-compressed DNG".into())), - 34892 => return Err(Error::Unsupported("lossy JPEG DNG".into())), + // DNG 1.4: lossy JPEG only for 8-bit LinearRaw (or mask) data + 34892 if photometric != PHOTOMETRIC_LINEAR_RAW => { + return Err(Error::Unsupported("lossy JPEG DNG of CFA data (DNG allows it only for 8-bit LinearRaw)".into())); + } + 34892 if bits != 8 => return Err(Error::Unsupported(format!("lossy JPEG DNG with {bits}-bit samples (DNG allows only 8)"))), + 34892 => {} 52546 => return Err(Error::Unsupported("JPEG XL DNG".into())), c => return Err(Error::Unsupported(format!("DNG compression {c}"))), } @@ -206,6 +214,7 @@ fn decode_chunk(c: &Chunk, l: &RawLayout, big_endian: bool) -> Result> let need = c.w * c.h * l.samples_per_pixel; match l.compression { 1 => unpack(c.bytes, c.w, c.h, l.samples_per_pixel, l.bits, big_endian), + 34892 => lossy_jpeg(c.bytes, need), _ => { // the JPEG's samples fill the chunk in raster order (its own width × components may // differ from the chunk's, e.g. two components of half the width) @@ -222,6 +231,38 @@ fn decode_chunk(c: &Chunk, l: &RawLayout, big_endian: bool) -> Result> } } +/// Decode a baseline (lossy) JPEG tile or strip (Compression 34892) into `need` 8-bit samples: +/// one component (greyscale) or three (colour, converted from YCbCr by the JPEG decoder as the +/// JPEG says), in raster order. +fn lossy_jpeg(bytes: &[u8], need: usize) -> Result> { + use image::ImageDecoder; + let bad = |_| Error::Malformed("lossy JPEG tile"); + let decode = || -> Result> { + let mut d = image::codecs::jpeg::JpegDecoder::new(std::io::Cursor::new(bytes)).map_err(bad)?; + let channels = match d.color_type() { + image::ColorType::L8 => 1usize, + image::ColorType::Rgb8 => 3, + c => return Err(Error::Unsupported(format!("lossy JPEG DNG tile of colour type {c:?}"))), + }; + let (w, h) = d.dimensions(); + let total = (w as usize).checked_mul(h as usize).and_then(|v| v.checked_mul(channels)).ok_or(Error::TooLarge)?; + if total < need { + return Err(Error::Malformed("lossy JPEG smaller than its tile or strip")); + } + if total > need.saturating_mul(4).saturating_add(65_536) || total > MAX_SAMPLES { + return Err(Error::TooLarge); + } + let mut limits = image::Limits::default(); + limits.max_alloc = Some(total as u64 * 4 + (1 << 20)); + d.set_limits(limits).map_err(bad)?; + let mut buf = vec![0u8; total]; + d.read_image(&mut buf).map_err(bad)?; + Ok(buf.iter().take(need).map(|&v| u16::from(v)).collect()) + }; + // a third-party decoder: a panic in it is an error for this file, not a crash + std::panic::catch_unwind(std::panic::AssertUnwindSafe(decode)).unwrap_or(Err(Error::Malformed("lossy JPEG decoder failed"))) +} + /// Read the raw image of `ifd`. pub fn read(t: &Tiff, ifd: &Ifd) -> Result { let l = layout(t, ifd)?; @@ -231,9 +272,12 @@ pub fn read(t: &Tiff, ifd: &Ifd) -> Result { return Err(Error::TooLarge); } let chunks = chunks(t, ifd, &l)?; - // every sample takes at least one bit of file (a Huffman code is at least one bit long) + // every sample takes at least one bit of file (a Huffman code is at least one bit long); + // lossy JPEG at least two bits per 8 × 8 block and component (a DC code and an end of + // block), so a tiny file cannot claim a huge picture + let per_byte = if l.compression == 34892 { 256 } else { 8 }; let coded: usize = chunks.iter().map(|c| c.bytes.len()).fold(0usize, usize::saturating_add); - if total.div_ceil(8) > coded { + if total.div_ceil(per_byte) > coded { return Err(Error::Truncated); } let mut decoded_samples = 0usize; @@ -243,23 +287,50 @@ pub fn read(t: &Tiff, ifd: &Ifd) -> Result { if n > MAX_SAMPLES || decoded_samples > MAX_DECODED_TILE_SAMPLES { return Err(Error::TooLarge); } - if n.div_ceil(8) > c.bytes.len() { + if n.div_ceil(per_byte) > c.bytes.len() { return Err(Error::Truncated); } } - let mut data = vec![0u16; total]; - let row_len = l.width * spp; - for c in &chunks { - // Copy each padded tile immediately. Retaining every decoded tile until assembly lets - // repeated ranges multiply memory beyond the final raw frame's bounded allocation. - let d = decode_chunk(c, &l, t.big_endian)?; - let cw = c.w.min(l.width - c.x); - let ch = c.h.min(l.height - c.y); - for r in 0..ch { - let src = d.get(r * c.w * spp..(r * c.w + cw) * spp).ok_or(Error::Truncated)?; - let at = (c.y + r) * row_len + c.x * spp; - data.get_mut(at..at + cw * spp).ok_or(Error::Truncated)?.copy_from_slice(src); + let raw = |data| RawImage { width: l.width, height: l.height, samples_per_pixel: spp, bits: l.bits, data }; + // one chunk covering the image (a single strip): its samples are the image, no copy + if let [c] = chunks.as_slice() + && (c.x, c.y, c.w) == (0, 0, l.width) + && c.h >= l.height + { + let mut d = decode_chunk(c, &l, t.big_endian)?; + d.truncate(total); + if d.len() == total { + return Ok(raw(d)); } + return Err(Error::Truncated); } - Ok(RawImage { width: l.width, height: l.height, samples_per_pixel: spp, bits: l.bits, data }) + // bands of one row of tiles (or one strip): decoded in parallel, each band copying its own + // chunks into place + let row_len = l.width * spp; + let (Some(first), true) = (chunks.first(), row_len > 0) else { return Err(Error::Malformed("no strips or tiles")) }; + let (band_rows, across) = (first.h.max(1), l.width.div_ceil(first.w.max(1))); + let mut data = vec![0u16; total]; + let band = |k: usize, out: &mut [u16]| -> Result<()> { + let mine = chunks.get(k * across..((k + 1) * across).min(chunks.len())).unwrap_or(&[]); + #[cfg(feature = "threads")] + let decoded: Vec>> = mine.par_iter().map(|c| decode_chunk(c, &l, t.big_endian)).collect(); + #[cfg(not(feature = "threads"))] + let decoded: Vec>> = mine.iter().map(|c| decode_chunk(c, &l, t.big_endian)).collect(); + for (c, d) in mine.iter().zip(decoded) { + let d = d?; + let cw = c.w.min(l.width.saturating_sub(c.x)); + let ch = c.h.min(l.height.saturating_sub(c.y)); + for r in 0..ch { + let src = d.get(r * c.w * spp..(r * c.w + cw) * spp).ok_or(Error::Truncated)?; + let at = r * row_len + c.x * spp; + out.get_mut(at..at + cw * spp).ok_or(Error::Truncated)?.copy_from_slice(src); + } + } + Ok(()) + }; + #[cfg(feature = "threads")] + data.par_chunks_mut(band_rows * row_len).enumerate().try_for_each(|(k, out)| band(k, out))?; + #[cfg(not(feature = "threads"))] + data.chunks_mut(band_rows * row_len).enumerate().try_for_each(|(k, out)| band(k, out))?; + Ok(raw(data)) } diff --git a/crates/dng/src/reference/demosaic.rs b/crates/dng/src/reference/demosaic.rs new file mode 100644 index 0000000000..3ef013b77a --- /dev/null +++ b/crates/dng/src/reference/demosaic.rs @@ -0,0 +1,131 @@ +//! The demosaic of PR #1068, kept unchanged as the reference (tests and benchmarks only). + +#[cfg(feature = "threads")] +use rayon::prelude::*; + +use crate::demosaic::Pattern; +use crate::options::DemosaicQuality; + +const PAD: usize = 2; + +/// Mirror index `i` (which may be up to PAD outside) into 0..n. +fn mirror(i: isize, n: usize) -> usize { + let n = n as isize; + let mut i = i; + if i < 0 { + i = -i; + } + if i >= n { + i = 2 * (n - 1) - i; + } + i.clamp(0, n - 1) as usize +} + +/// The plane extended by PAD mirrored pixels on each side. +fn padded(plane: &[f32], w: usize, h: usize) -> (Vec, usize) { + let pw = w + 2 * PAD; + let ph = h + 2 * PAD; + let mut out = vec![0f32; pw * ph]; + for y in 0..ph { + let sy = mirror(y as isize - PAD as isize, h); + for x in 0..pw { + let sx = mirror(x as isize - PAD as isize, w); + out[y * pw + x] = plane.get(sy * w + sx).copied().unwrap_or(0.0); + } + } + (out, pw) +} + +/// Demosaic `plane` (`w` × `h`, one sample per pixel) into interleaved RGB. +pub fn demosaic(plane: &[f32], w: usize, h: usize, pattern: &Pattern) -> Vec { + demosaic_with(plane, w, h, pattern, DemosaicQuality::High) +} + +/// [`demosaic`] at a chosen quality: [`DemosaicQuality::Fast`] is bilinear for every pattern. +pub fn demosaic_with(plane: &[f32], w: usize, h: usize, pattern: &Pattern, quality: DemosaicQuality) -> Vec { + let mut out = vec![0f32; w * h * 3]; + if w == 0 || h == 0 || plane.len() < w * h { + return out; + } + let (p, pw) = padded(plane, w, h); + let bayer = quality == DemosaicQuality::High && pattern.is_bayer() && w >= 3 && h >= 3; + let row = |y: usize, o: &mut [f32]| { + if bayer { + malvar_row(&p, pw, y, w, pattern, o); + } else { + bilinear_row(&p, pw, y, w, h, pattern, o); + } + }; + #[cfg(feature = "threads")] + out.par_chunks_mut(w * 3).enumerate().for_each(|(y, o)| row(y, o)); + #[cfg(not(feature = "threads"))] + out.chunks_mut(w * 3).enumerate().for_each(|(y, o)| row(y, o)); + out +} + +fn malvar_row(p: &[f32], pw: usize, y: usize, w: usize, pat: &Pattern, out: &mut [f32]) { + let py = y + PAD; + for x in 0..w { + let px = x + PAD; + let at = |dx: isize, dy: isize| p[((py as isize + dy) as usize) * pw + (px as isize + dx) as usize]; + let c = at(0, 0); + let cross1 = at(-1, 0) + at(1, 0) + at(0, -1) + at(0, 1); + let cross2 = at(-2, 0) + at(2, 0) + at(0, -2) + at(0, 2); + let diag = at(-1, -1) + at(1, -1) + at(-1, 1) + at(1, 1); + let me = pat.color(x, y); + let rgb = match me { + 1 => { + let h1 = at(-1, 0) + at(1, 0); + let h2 = at(-2, 0) + at(2, 0); + let v1 = at(0, -1) + at(0, 1); + let v2 = at(0, -2) + at(0, 2); + // colour along the row and along the column + let along_row = (5.0 * c + 4.0 * h1 - h2 - diag + 0.5 * v2) / 8.0; + let along_col = (5.0 * c + 4.0 * v1 - v2 - diag + 0.5 * h2) / 8.0; + if pat.color(x + 1, y) == 0 { [along_row, c, along_col] } else { [along_col, c, along_row] } + } + _ => { + let g = (4.0 * c + 2.0 * cross1 - cross2) / 8.0; + let other = (6.0 * c + 2.0 * diag - 1.5 * cross2) / 8.0; + if me == 0 { [c, g, other] } else { [other, g, c] } + } + }; + if let Some(o) = out.get_mut(x * 3..x * 3 + 3) { + o.copy_from_slice(&rgb); + } + } +} + +fn bilinear_row(p: &[f32], pw: usize, y: usize, w: usize, h: usize, pat: &Pattern, out: &mut [f32]) { + for x in 0..w { + let me = pat.color(x, y); + let mut rgb = [0f32; 3]; + for (k, v) in rgb.iter_mut().enumerate() { + if usize::from(me) == k { + *v = p[(y + PAD) * pw + x + PAD]; + continue; + } + for r in 1..=PAD as isize { + let (mut sum, mut n) = (0f32, 0u32); + for dy in -r..=r { + for dx in -r..=r { + // the colour of a mirrored neighbour is the colour where it was mirrored from + let sx = mirror(x as isize + dx, w); + let sy = mirror(y as isize + dy, h); + if usize::from(pat.color(sx, sy)) == k { + sum += p[(sy + PAD) * pw + sx + PAD]; + n += 1; + } + } + } + if n > 0 { + *v = sum / n as f32; + break; + } + } + } + if let Some(o) = out.get_mut(x * 3..x * 3 + 3) { + o.copy_from_slice(&rgb); + } + } +} diff --git a/crates/dng/src/reference/ljpeg.rs b/crates/dng/src/reference/ljpeg.rs new file mode 100644 index 0000000000..8b3fa935d3 --- /dev/null +++ b/crates/dng/src/reference/ljpeg.rs @@ -0,0 +1,384 @@ +//! The lossless JPEG decoder of PR #1068, kept unchanged as the reference (tests and benchmarks +//! only). + +use crate::{Error, Result}; + +/// A decoded lossless JPEG image. +#[derive(Debug)] +pub struct Decoded { + pub width: usize, + pub height: usize, + pub components: usize, + pub precision: u8, + pub predictor: u8, + /// Samples, interleaved by component, line after line, shifted back up by the point transform. + pub samples: Vec, +} + +const LOOKUP_BITS: u32 = 9; + +/// A Huffman table for decoding (T.81 Annex C and F.2.2.3). +#[derive(Clone)] +struct Huffman { + /// For codes of up to LOOKUP_BITS bits: (length, value) indexed by the next LOOKUP_BITS bits; + /// length 0 means "longer code". + fast: Vec<(u8, u8)>, + maxcode: [i32; 17], + valptr: [i32; 17], + mincode: [i32; 17], + values: Vec, +} + +impl Huffman { + fn new(bits: &[u8; 16], values: Vec) -> Result { + // C.2: code sizes and codes + let mut sizes = Vec::new(); + for (l, &n) in bits.iter().enumerate() { + for _ in 0..n { + sizes.push(l as u8 + 1); + } + } + if sizes.len() != values.len() || sizes.is_empty() { + return Err(Error::Malformed("lossless JPEG: Huffman table counts")); + } + let mut codes = Vec::with_capacity(sizes.len()); + let mut code: u32 = 0; + let mut si = sizes[0]; + let mut k = 0; + while k < sizes.len() { + while k < sizes.len() && sizes[k] == si { + codes.push(code); + code += 1; + k += 1; + } + // the codes of one length must fit in that length + if code > (1u32 << si) { + return Err(Error::Malformed("lossless JPEG: Huffman table over-subscribed")); + } + code <<= 1; + si += 1; + } + // F.2.2.3 decoder tables + let mut maxcode = [-1i32; 17]; + let mut valptr = [0i32; 17]; + let mut mincode = [0i32; 17]; + let mut j = 0usize; + for l in 1..=16usize { + let n = usize::from(bits[l - 1]); + if n == 0 { + maxcode[l] = -1; + } else { + valptr[l] = j as i32; + mincode[l] = codes[j] as i32; + j += n; + maxcode[l] = codes[j - 1] as i32; + } + } + let mut fast = vec![(0u8, 0u8); 1 << LOOKUP_BITS]; + for (i, (&c, &s)) in codes.iter().zip(&sizes).enumerate() { + let s = u32::from(s); + if s <= LOOKUP_BITS { + let shift = LOOKUP_BITS - s; + let base = (c << shift) as usize; + for f in fast.iter_mut().skip(base).take(1 << shift) { + *f = (s as u8, values[i]); + } + } + } + Ok(Huffman { fast, maxcode, valptr, mincode, values }) + } +} + +/// Entropy-coded data reader: removes stuffed zero bytes and stops at markers (feeding zeros). +struct Bits<'a> { + data: &'a [u8], + pos: usize, + acc: u64, + n: u32, + /// The marker that ended the data (after 0xFF), if one was met. + marker: Option, +} + +impl<'a> Bits<'a> { + fn new(data: &'a [u8], pos: usize) -> Self { + Bits { data, pos, acc: 0, n: 0, marker: None } + } + + fn fill(&mut self) { + while self.n <= 56 { + let mut byte = 0u8; + if self.marker.is_none() { + match self.data.get(self.pos) { + Some(&0xFF) => match self.data.get(self.pos + 1) { + Some(&0x00) => { + byte = 0xFF; + self.pos += 2; + } + Some(&m) => self.marker = Some(m), + None => self.marker = Some(0xD9), + }, + Some(&b) => { + byte = b; + self.pos += 1; + } + None => self.marker = Some(0xD9), + } + } + self.acc |= u64::from(byte) << (56 - self.n); + self.n += 8; + } + } + + fn peek(&mut self, k: u32) -> u32 { + if self.n < k { + self.fill(); + } + (self.acc >> (64 - k)) as u32 + } + + fn skip(&mut self, k: u32) { + self.acc <<= k; + self.n -= k; + } + + fn bits(&mut self, k: u32) -> u32 { + if k == 0 { + return 0; + } + let v = self.peek(k); + self.skip(k); + v + } + + /// Decode one Huffman symbol (F.2.2.3). + fn decode(&mut self, h: &Huffman) -> Result { + let p = self.peek(LOOKUP_BITS) as usize; + if let Some(&(len, v)) = h.fast.get(p) + && len > 0 + { + self.skip(u32::from(len)); + return Ok(v); + } + let mut code = self.bits(1) as i32; + for l in 1..=16usize { + if code <= h.maxcode[l] { + let i = h.valptr[l] + code - h.mincode[l]; + return h.values.get(i as usize).copied().ok_or(Error::Malformed("lossless JPEG: bad Huffman code")); + } + code = (code << 1) | self.bits(1) as i32; + } + Err(Error::Malformed("lossless JPEG: bad Huffman code")) + } + + /// Skip to the byte after the next RSTn marker, starting a new bit stream there. + fn restart(&mut self) -> Result<()> { + // the decoder may have looked ahead into fill bits; find the marker in the bytes + let mut p = self.pos; + loop { + match (self.data.get(p), self.data.get(p + 1)) { + (Some(0xFF), Some(&m)) if (0xD0..=0xD7).contains(&m) => { + *self = Bits::new(self.data, p + 2); + return Ok(()); + } + (Some(_), Some(_)) => p += 1, + _ => return Err(Error::Malformed("lossless JPEG: missing restart marker")), + } + } + } +} + +struct Component { + id: u8, + table: usize, +} + +/// Decode a lossless JPEG stream. `max_samples` caps the picture (width × height × components). +pub fn decode(data: &[u8], max_samples: usize) -> Result { + if data.get(0..2) != Some(&[0xFF, 0xD8]) { + return Err(Error::Malformed("lossless JPEG: no SOI marker")); + } + let mut pos = 2usize; + let mut tables: [Option; 4] = [None, None, None, None]; + let mut frame: Option<(u8, usize, usize, Vec)> = None; + let mut restart_interval = 0usize; + loop { + // markers may be preceded by fill bytes 0xFF + while data.get(pos) == Some(&0xFF) && data.get(pos + 1) == Some(&0xFF) { + pos += 1; + } + let (Some(&0xFF), Some(&m)) = (data.get(pos), data.get(pos + 1)) else { + return Err(Error::Malformed("lossless JPEG: expected a marker")); + }; + pos += 2; + if m == 0xD9 { + return Err(Error::Malformed("lossless JPEG: no scan")); + } + if (0xD0..=0xD8).contains(&m) || m == 0x01 { + continue; // markers without a segment + } + let len_b = data.get(pos..pos + 2).ok_or(Error::Truncated)?; + let len = usize::from(u16::from_be_bytes([len_b[0], len_b[1]])); + if len < 2 { + return Err(Error::Malformed("lossless JPEG: segment length")); + } + let seg = data.get(pos + 2..pos + len).ok_or(Error::Truncated)?; + match m { + 0xC3 => { + let [p, y0, y1, x0, x1, nf, rest @ ..] = seg else { + return Err(Error::Malformed("lossless JPEG: short frame header")); + }; + let (precision, nf) = (*p, usize::from(*nf)); + let height = usize::from(u16::from_be_bytes([*y0, *y1])); + let width = usize::from(u16::from_be_bytes([*x0, *x1])); + if !(2..=16).contains(&precision) { + return Err(Error::Unsupported(format!("lossless JPEG precision {precision}"))); + } + if height == 0 { + return Err(Error::Unsupported("lossless JPEG with the height in a DNL marker".into())); + } + if width == 0 || !(1..=4).contains(&nf) || rest.len() < nf * 3 { + return Err(Error::Malformed("lossless JPEG: frame header")); + } + let mut comps = Vec::with_capacity(nf); + for c in rest.as_chunks::<3>().0.iter().take(nf) { + if c[1] != 0x11 { + return Err(Error::Unsupported("subsampled lossless JPEG components".into())); + } + comps.push(Component { id: c[0], table: 0 }); + } + frame = Some((precision, width, height, comps)); + } + 0xC4 => { + let mut s = seg; + while !s.is_empty() { + let [tcth, counts @ ..] = s else { break }; + let bits: [u8; 16] = counts.get(..16).and_then(|b| b.try_into().ok()).ok_or(Error::Truncated)?; + let n: usize = bits.iter().map(|&b| usize::from(b)).sum(); + let values = counts.get(16..16 + n).ok_or(Error::Truncated)?.to_vec(); + let th = usize::from(tcth & 0x0F); + if tcth >> 4 != 0 || th > 3 { + return Err(Error::Malformed("lossless JPEG: Huffman table class or id")); + } + if values.iter().any(|&v| v > 16) { + return Err(Error::Malformed("lossless JPEG: difference category above 16")); + } + tables[th] = Some(Huffman::new(&bits, values)?); + s = counts.get(16 + n..).unwrap_or(&[]); + } + } + 0xDD => { + let r = seg.get(0..2).ok_or(Error::Truncated)?; + restart_interval = usize::from(u16::from_be_bytes([r[0], r[1]])); + } + 0xDA => { + let Some((precision, width, height, mut comps)) = frame.take() else { + return Err(Error::Malformed("lossless JPEG: scan before the frame header")); + }; + let ns = usize::from(*seg.first().ok_or(Error::Truncated)?); + if ns != comps.len() { + return Err(Error::Unsupported("lossless JPEG with more than one scan".into())); + } + for k in 0..ns { + let c = seg.get(1 + 2 * k..3 + 2 * k).ok_or(Error::Truncated)?; + let comp = comps.iter_mut().find(|x| x.id == c[0]).ok_or(Error::Malformed("lossless JPEG: unknown scan component"))?; + comp.table = usize::from(c[1] >> 4); + if comp.table > 3 { + return Err(Error::Malformed("lossless JPEG: Huffman table id")); + } + } + let tail = seg.get(1 + 2 * ns..4 + 2 * ns).ok_or(Error::Truncated)?; + let (predictor, pt) = (tail[0], tail[2] & 0x0F); + if !(1..=7).contains(&predictor) { + return Err(Error::Unsupported(format!("lossless JPEG predictor {predictor}"))); + } + if pt >= precision { + return Err(Error::Malformed("lossless JPEG: point transform")); + } + let n = width.checked_mul(height).and_then(|v| v.checked_mul(comps.len())).ok_or(Error::TooLarge)?; + if n > max_samples { + return Err(Error::TooLarge); + } + let mut hufs = Vec::with_capacity(comps.len()); + for c in &comps { + hufs.push(tables.get(c.table).and_then(|t| t.as_ref()).ok_or(Error::Malformed("lossless JPEG: missing Huffman table"))?); + } + let samples = scan(data, pos + len, &hufs, width, height, precision, predictor, pt, restart_interval)?; + return Ok(Decoded { width, height, components: comps.len(), precision, predictor, samples }); + } + 0xC0..=0xC2 | 0xC5..=0xC7 | 0xC9..=0xCB | 0xCD..=0xCF => { + return Err(Error::Unsupported("lossy (DCT) JPEG in a DNG".into())); + } + _ => {} // APPn, COM, DQT … are skipped + } + pos += len; + } +} + +/// Decode the entropy-coded data of one scan (H.1.2). +fn scan(data: &[u8], start: usize, hufs: &[&Huffman], width: usize, height: usize, precision: u8, predictor: u8, pt: u8, restart: usize) -> Result> { + let nc = hufs.len(); + let row_len = width * nc; + if restart != 0 && !restart.is_multiple_of(width) { + return Err(Error::Unsupported("lossless JPEG restart interval that splits a line".into())); + } + let lines_per_interval = if restart == 0 { usize::MAX } else { restart / width }; + let mut out = vec![0u16; row_len * height]; + let mut bits = Bits::new(data, start); + let initial = 1i32 << (precision - pt - 1); + let mut line_in_interval = 0usize; + for y in 0..height { + if line_in_interval == lines_per_interval { + bits.restart()?; + line_in_interval = 0; + } + let first_line = line_in_interval == 0; + let (above, cur) = if y == 0 { + (None, &mut out[..row_len]) + } else { + let (a, c) = out.split_at_mut(y * row_len); + (Some(&a[(y - 1) * row_len..]), &mut c[..row_len]) + }; + for x in 0..width { + for (c, h) in hufs.iter().enumerate() { + let i = x * nc + c; + let ssss = u32::from(bits.decode(h)?); + let diff = match ssss { + 0 => 0, + 16 => 32768, + s => { + let v = bits.bits(s) as i32; + if v < (1 << (s - 1)) { v - (1 << s) + 1 } else { v } + } + }; + let px = match above { + _ if first_line && x == 0 => initial, + _ if first_line => i32::from(cur[i - nc]), + Some(a) if x == 0 => i32::from(a[i]), + Some(a) => predict(predictor, i32::from(cur[i - nc]), i32::from(a[i]), i32::from(a[i - nc])), + None => initial, + }; + cur[i] = ((px + diff) & 0xFFFF) as u16; + } + } + line_in_interval += 1; + } + if pt > 0 { + for s in &mut out { + *s <<= pt; + } + } + Ok(out) +} + +/// Table H.1 predictors (`ra` left, `rb` above, `rc` above left). +pub fn predict(p: u8, ra: i32, rb: i32, rc: i32) -> i32 { + match p { + 1 => ra, + 2 => rb, + 3 => rc, + 4 => ra + rb - rc, + 5 => ra + ((rb - rc) >> 1), + 6 => rb + ((ra - rc) >> 1), + _ => (ra + rb) >> 1, + } +} diff --git a/crates/dng/src/reference/mod.rs b/crates/dng/src/reference/mod.rs new file mode 100644 index 0000000000..9236526b5c --- /dev/null +++ b/crates/dng/src/reference/mod.rs @@ -0,0 +1,44 @@ +//! The decoder and pipeline of PR #1068, kept unchanged (feature `writer`, tests and benchmarks +//! only) so the faster implementation can be checked for bit-identical output and benchmarked +//! against it in the same run. + +#![allow(clippy::too_many_arguments)] + +pub mod demosaic; +pub mod ljpeg; +pub mod process; +pub mod raw; + +use crate::options::{ProcessOptions, WhiteBalance}; +use crate::{Dng, Error, Image, Result, color}; + +impl Dng<'_> { + /// The raw samples as PR #1068 read them. + pub fn reference_raw(&self) -> Result { + raw::read(&self.tiff, self.raw_ifd()) + } + + /// [`Dng::develop_with`] as PR #1068 (plus the `ProcessOptions` commit) developed. + pub fn reference_develop(&self, o: &ProcessOptions) -> Result { + let info = self.info()?; + let p = self.params()?; + let raw = self.reference_raw()?; + let mut calibration = self.calibration()?; + calibration.baseline_exposure += o.exposure_stops(); + if let WhiteBalance::Custom { temperature, tint } = o.white_balance + && calibration.color_matrix.iter().any(Option::is_some) + { + let xy = color::temperature_tint_to_xy(temperature, tint).ok_or(Error::Unsupported("a non-finite white balance".into()))?; + calibration.as_shot_neutral = None; + calibration.as_shot_white_xy = Some(xy); + } + let solved = if p.cfa.is_none() && raw.samples_per_pixel == 1 { + let e = 2f64.powf(calibration.baseline_exposure.clamp(-16.0, 16.0)); + color::Solved { neutral: [1.0; 3], matrix: color::scale(&color::IDENTITY, e), white_xy: color::D65 } + } else { + calibration.solve()? + }; + let d = process::develop_with(&raw, &p, &solved.neutral, &solved.matrix, o)?; + Ok(Image { info, width: d.width, height: d.height, rgba: d.rgba }) + } +} diff --git a/crates/dng/src/reference/process.rs b/crates/dng/src/reference/process.rs new file mode 100644 index 0000000000..13018717fa --- /dev/null +++ b/crates/dng/src/reference/process.rs @@ -0,0 +1,118 @@ +//! The processing of PR #1068 (normalise, demosaic, crop, matrix), kept unchanged as the reference +//! the faster pipeline is compared against (tests and benchmarks only). + +#[cfg(feature = "threads")] +use rayon::prelude::*; + +use super::demosaic::demosaic_with; +use crate::color::{M3, V3}; +use crate::options::ProcessOptions; +use crate::process::{Developed, Params, active_rect, crop_rect}; + +type Rect = (usize, usize, usize, usize); +use crate::raw::RawImage; +use crate::{Error, Result}; + +/// Normalised, white-balanced samples of the active area (steps 1–4), one plane per sample. +fn normalise(raw: &RawImage, p: &Params, (ax, ay, aw, ah): Rect, neutral: &V3, planes: usize) -> Result> { + let spp = raw.samples_per_pixel; + let (br, bc) = (p.black_repeat.0.max(1), p.black_repeat.1.max(1)); + if br > 64 || bc > 64 { + return Err(Error::Unsupported("a BlackLevelRepeatDim larger than 64".into())); + } + let need = br * bc * spp; + let black: Vec = match p.black.len() { + 0 => vec![0.0; need], + 1 => vec![p.black[0]; need], + n if n >= need => p.black[..need].to_vec(), + _ => return Err(Error::Malformed("BlackLevel has too few values")), + }; + let max_code = ((1u32 << raw.bits) - 1) as f64; + let white: Vec = (0..spp).map(|s| p.white.get(s).or(p.white.first()).copied().unwrap_or(max_code)).collect(); + let lin = p.linearization.as_deref().filter(|t| !t.is_empty()); + let inv_neutral: Vec = (0..planes).map(|c| (1.0 / neutral.get(c).copied().unwrap_or(1.0)) as f32).collect(); + let cfa = p.cfa.as_ref(); + let mut out = vec![0f32; aw * ah * spp]; + let row = |y: usize, o: &mut [f32]| { + let dv = p.black_delta_v.get(y).copied().unwrap_or(0.0); + let src_row = (ay + y) * raw.width * spp; + for x in 0..aw { + let dh = p.black_delta_h.get(x).copied().unwrap_or(0.0); + for s in 0..spp { + let code = raw.data.get(src_row + (ax + x) * spp + s).copied().unwrap_or(0); + let v = match lin { + Some(t) => f64::from(t.get(usize::from(code)).or(t.last()).copied().unwrap_or(0)), + None => f64::from(code), + }; + let b = black[((y % br) * bc + x % bc) * spp + s] + dh + dv; + let range = white[s] - b; + let n = if range > 0.0 { ((v - b) / range).clamp(0.0, 1.0) } else { 0.0 }; + let plane = match cfa { + Some(c) => usize::from(c.color(x, y)), + None => s, + }; + let wb = inv_neutral.get(plane).copied().unwrap_or(1.0); + if let Some(o) = o.get_mut(x * spp + s) { + *o = (n as f32 * wb).min(1.0); + } + } + } + }; + #[cfg(feature = "threads")] + out.par_chunks_mut((aw * spp).max(1)).enumerate().for_each(|(y, o)| row(y, o)); + #[cfg(not(feature = "threads"))] + out.chunks_mut((aw * spp).max(1)).enumerate().for_each(|(y, o)| row(y, o)); + Ok(out) +} + +/// Develop `raw` (steps 1–7). `neutral` and `matrix` come from [`crate::color::Calibration::solve`]. +pub fn develop(raw: &RawImage, p: &Params, neutral: &V3, matrix: &M3) -> Result { + develop_with(raw, p, neutral, matrix, &ProcessOptions::default()) +} + +/// [`develop`] with raw settings. White balance and exposure are already in `neutral` and +/// `matrix` ([`crate::Dng::develop_with`] solves them); this applies the demosaic quality. +pub fn develop_with(raw: &RawImage, p: &Params, neutral: &V3, matrix: &M3, o: &ProcessOptions) -> Result { + let spp = raw.samples_per_pixel; + if p.cfa.is_some() && spp != 1 { + return Err(Error::Malformed("a CFA image with more than one sample per pixel")); + } + if p.cfa.is_none() && !matches!(spp, 1 | 3) { + return Err(Error::Unsupported(format!("LinearRaw with {spp} samples per pixel"))); + } + let active = active_rect(p, raw.width, raw.height)?; + let (_, _, aw, ah) = active; + let planes = if p.cfa.is_some() { 3 } else { spp }; + let norm = normalise(raw, p, active, neutral, planes)?; + // three samples per pixel + let rgb = match (&p.cfa, spp) { + (Some(pat), _) => demosaic_with(&norm, aw, ah, pat, o.demosaic), + (None, 3) => norm, + _ => norm.iter().flat_map(|&v| [v, v, v]).collect(), + }; + let (cx, cy, cw, ch) = crop_rect(p, aw, ah); + let m: [[f32; 3]; 3] = matrix.map(|r| r.map(|v| v as f32)); + let mono = p.cfa.is_none() && spp == 1; + let mut rgba = vec![0f32; cw * ch * 4]; + let row = |y: usize, o: &mut [f32]| { + let src = (cy + y) * aw + cx; + for x in 0..cw { + let Some(c) = rgb.get((src + x) * 3..(src + x) * 3 + 3) else { continue }; + let px = if mono { + // grey: the matrix's luminance gain only + let g = c[0] * (m[1][0] + m[1][1] + m[1][2]); + [g, g, g] + } else { + [0, 1, 2].map(|i| m[i][0] * c[0] + m[i][1] * c[1] + m[i][2] * c[2]) + }; + if let Some(d) = o.get_mut(x * 4..x * 4 + 4) { + d.copy_from_slice(&[px[0].max(0.0), px[1].max(0.0), px[2].max(0.0), 1.0]); + } + } + }; + #[cfg(feature = "threads")] + rgba.par_chunks_mut(cw * 4).enumerate().for_each(|(y, o)| row(y, o)); + #[cfg(not(feature = "threads"))] + rgba.chunks_mut(cw * 4).enumerate().for_each(|(y, o)| row(y, o)); + Ok(Developed { width: cw, height: ch, rgba }) +} diff --git a/crates/dng/src/reference/raw.rs b/crates/dng/src/reference/raw.rs new file mode 100644 index 0000000000..0b14c6992a --- /dev/null +++ b/crates/dng/src/reference/raw.rs @@ -0,0 +1,156 @@ +//! The raw sample reading of PR #1068, kept unchanged as the reference the faster reader is +//! compared against (tests and benchmarks only). + +#[cfg(feature = "threads")] +use rayon::prelude::*; + +use super::ljpeg; +use crate::raw::{MAX_SAMPLES, RawImage, RawLayout, layout}; +use crate::tiff::{Ifd, Tiff, tags}; +use crate::{Error, Result}; + +fn uint_or(t: &Tiff, ifd: &Ifd, tag: u16, default: u64) -> Result { + match ifd.get(tag) { + Some(e) => t.uint(e), + None => Ok(default), + } +} + +/// One strip or tile: where it sits in the image and its bytes. +struct Chunk<'a> { + x: usize, + y: usize, + /// Pixel size of the chunk as coded (tiles are whole even at the image edge). + w: usize, + h: usize, + bytes: &'a [u8], +} + +fn chunks<'a>(t: &Tiff<'a>, ifd: &Ifd, l: &RawLayout) -> Result>> { + let (offs, counts, tw, th) = if let (Some(o), Some(c)) = (ifd.get(tags::TILE_OFFSETS), ifd.get(tags::TILE_BYTE_COUNTS)) { + let tw = uint_or(t, ifd, tags::TILE_WIDTH, 0)? as usize; + let th = uint_or(t, ifd, tags::TILE_LENGTH, 0)? as usize; + if tw == 0 || th == 0 || tw > 65_535 || th > 65_535 { + return Err(Error::Malformed("tile size")); + } + (t.uints(o)?, t.uints(c)?, tw, th) + } else if let (Some(o), Some(c)) = (ifd.get(tags::STRIP_OFFSETS), ifd.get(tags::STRIP_BYTE_COUNTS)) { + let rows = (uint_or(t, ifd, tags::ROWS_PER_STRIP, l.height as u64)? as usize).clamp(1, l.height); + (t.uints(o)?, t.uints(c)?, l.width, rows) + } else { + return Err(Error::Malformed("no strips or tiles")); + }; + let across = l.width.div_ceil(tw); + let down = l.height.div_ceil(th); + let n = across * down; + if offs.len() < n || counts.len() < n { + return Err(Error::Malformed("fewer strips or tiles than the image needs")); + } + let mut out = Vec::with_capacity(n); + for i in 0..n { + let (x, y) = ((i % across) * tw, (i / across) * th); + let start = usize::try_from(offs[i]).map_err(|_| Error::Truncated)?; + let len = usize::try_from(counts[i]).map_err(|_| Error::Truncated)?; + let end = start.checked_add(len).ok_or(Error::Truncated)?; + let bytes = t.data.get(start..end).ok_or(Error::Truncated)?; + // strips end at the image; tiles are coded whole + let h = if ifd.get(tags::TILE_OFFSETS).is_some() { th } else { th.min(l.height - y) }; + out.push(Chunk { x, y, w: tw, h, bytes }); + } + Ok(out) +} + +/// Unpack `w` × `h` × `spp` uncompressed samples of `bits` bits: 8 and 16-bit samples as bytes / +/// words in the file's byte order, other depths packed most significant bit first with each row +/// starting on a byte boundary (TIFF 6.0 §7, DNG BitsPerSample). +pub fn unpack(bytes: &[u8], w: usize, h: usize, spp: usize, bits: u32, big_endian: bool) -> Result> { + let per_row = w * spp; + let row_bytes = (per_row * bits as usize).div_ceil(8); + let need = row_bytes.checked_mul(h).ok_or(Error::TooLarge)?; + if bytes.len() < need { + return Err(Error::Truncated); + } + let mut out = Vec::with_capacity(per_row * h); + for row in bytes.chunks_exact(row_bytes.max(1)).take(h) { + match bits { + 8 => out.extend(row.iter().map(|&b| u16::from(b))), + 16 => out.extend(row.as_chunks::<2>().0.iter().map(|&b| if big_endian { u16::from_be_bytes(b) } else { u16::from_le_bytes(b) })), + _ => { + let mut acc: u32 = 0; + let mut n: u32 = 0; + let mut it = row.iter(); + for _ in 0..per_row { + while n < bits { + acc = (acc << 8) | u32::from(*it.next().ok_or(Error::Truncated)?); + n += 8; + } + n -= bits; + out.push(((acc >> n) & ((1 << bits) - 1)) as u16); + acc &= (1 << n) - 1; + } + } + } + } + Ok(out) +} + +/// Decode one chunk into `w` × `h` × `spp` samples. +fn decode_chunk(c: &Chunk, l: &RawLayout, big_endian: bool) -> Result> { + let need = c.w * c.h * l.samples_per_pixel; + match l.compression { + 1 => unpack(c.bytes, c.w, c.h, l.samples_per_pixel, l.bits, big_endian), + _ => { + // the JPEG's samples fill the chunk in raster order (its own width × components may + // differ from the chunk's, e.g. two components of half the width) + let d = ljpeg::decode(c.bytes, need.saturating_add(65_536 * 4))?; + if d.samples.len() < need { + return Err(Error::Malformed("lossless JPEG smaller than its tile or strip")); + } + let mut s = d.samples; + s.truncate(need); + Ok(s) + } + } +} + +/// Read the raw image of `ifd`. +pub fn read(t: &Tiff, ifd: &Ifd) -> Result { + let l = layout(t, ifd)?; + let spp = l.samples_per_pixel; + let total = l.width.checked_mul(l.height).and_then(|v| v.checked_mul(spp)).ok_or(Error::TooLarge)?; + if total > MAX_SAMPLES { + return Err(Error::TooLarge); + } + let chunks = chunks(t, ifd, &l)?; + // every sample takes at least one bit of file (a Huffman code is at least one bit long) + let coded: usize = chunks.iter().map(|c| c.bytes.len()).fold(0usize, usize::saturating_add); + if total / 8 > coded { + return Err(Error::Truncated); + } + for c in &chunks { + let n = c.w.checked_mul(c.h).and_then(|v| v.checked_mul(spp)).ok_or(Error::TooLarge)?; + if n > MAX_SAMPLES { + return Err(Error::TooLarge); + } + if n / 8 > c.bytes.len() { + return Err(Error::Truncated); + } + } + #[cfg(feature = "threads")] + let decoded: Vec>> = chunks.par_iter().map(|c| decode_chunk(c, &l, t.big_endian)).collect(); + #[cfg(not(feature = "threads"))] + let decoded: Vec>> = chunks.iter().map(|c| decode_chunk(c, &l, t.big_endian)).collect(); + let mut data = vec![0u16; total]; + let row_len = l.width * spp; + for (c, d) in chunks.iter().zip(decoded) { + let d = d?; + let cw = c.w.min(l.width - c.x); + let ch = c.h.min(l.height - c.y); + for r in 0..ch { + let src = d.get(r * c.w * spp..(r * c.w + cw) * spp).ok_or(Error::Truncated)?; + let at = (c.y + r) * row_len + c.x * spp; + data.get_mut(at..at + cw * spp).ok_or(Error::Truncated)?.copy_from_slice(src); + } + } + Ok(RawImage { width: l.width, height: l.height, samples_per_pixel: spp, bits: l.bits, data }) +} diff --git a/crates/dng/src/tiff.rs b/crates/dng/src/tiff.rs index 1e8e79c425..8b35df147e 100644 --- a/crates/dng/src/tiff.rs +++ b/crates/dng/src/tiff.rs @@ -331,6 +331,9 @@ pub mod tags { pub const CALIBRATION_ILLUMINANT_1: u16 = 50778; pub const CALIBRATION_ILLUMINANT_2: u16 = 50779; pub const ACTIVE_AREA: u16 = 50829; + pub const OPCODE_LIST_1: u16 = 51008; + pub const OPCODE_LIST_2: u16 = 51009; + pub const OPCODE_LIST_3: u16 = 51022; pub const FORWARD_MATRIX_1: u16 = 50964; pub const FORWARD_MATRIX_2: u16 = 50965; /// CinemaDNG: SMPTE ST 12-1 time codes, 8 bytes each. diff --git a/crates/dng/src/writer.rs b/crates/dng/src/writer.rs index 24203b1615..1bbb67b3d0 100644 --- a/crates/dng/src/writer.rs +++ b/crates/dng/src/writer.rs @@ -79,6 +79,11 @@ pub enum Compression { restart_lines: usize, pt: u8, }, + /// Baseline (lossy) JPEG at `quality` 1–100 (8-bit samples; one sample per pixel as a + /// greyscale JPEG, three as a colour one, which the encoder stores as YCbCr). + LossyJpeg { + quality: u8, + }, } #[derive(Clone, Copy, Debug, PartialEq, Eq)] @@ -318,6 +323,20 @@ fn write_entropy(out: Vec, coded: &[Coded], tables: &[HuffTable], restart_li out } +/// Encode `w` × `h` × `spp` 8-bit samples (1 or 3 per pixel) as a baseline JPEG (empty if the +/// encoder refuses them). +pub fn lossy_jpeg(samples: &[u16], w: usize, h: usize, spp: usize, quality: u8) -> Vec { + use image::ImageEncoder; + let bytes: Vec = samples.iter().map(|&v| v.min(255) as u8).collect(); + let color = if spp == 3 { image::ExtendedColorType::Rgb8 } else { image::ExtendedColorType::L8 }; + let mut out = Vec::new(); + let enc = image::codecs::jpeg::JpegEncoder::new_with_quality(&mut out, quality.clamp(1, 100)); + if enc.write_image(&bytes, w as u32, h as u32, color).is_err() { + return Vec::new(); + } + out +} + /// Encode `samples` (`width` × `height` × `comps`, interleaved) as lossless JPEG. pub fn lj92(samples: &[u16], params: &Lj92Params) -> Vec { let p = Lj92Params { comps: params.comps.clamp(1, 4), ..*params }; @@ -401,6 +420,7 @@ fn put_chunks(buf: &mut Vec, samples: &[u16], o: &Options) -> (Vec, Vec let params = Lj92Params { width: cw * spp / comps, height: rows, comps, precision: o.bits as u8, predictor, restart_lines, pt }; lj92(&chunk, ¶ms) } + Compression::LossyJpeg { quality } => lossy_jpeg(&chunk, cw, rows, spp, quality), }; if buf.len() % 2 == 1 { buf.push(0); @@ -423,7 +443,14 @@ pub fn write(samples: &[u16], o: &Options) -> Vec { (tags::IMAGE_WIDTH, Value::Long(vec![o.width as u32])), (tags::IMAGE_LENGTH, Value::Long(vec![o.height as u32])), (tags::BITS_PER_SAMPLE, Value::Short(vec![o.bits as u16; o.samples_per_pixel])), - (tags::COMPRESSION, Value::Short(vec![if o.compression == Compression::None { 1 } else { 7 }])), + ( + tags::COMPRESSION, + Value::Short(vec![match o.compression { + Compression::None => 1, + Compression::LosslessJpeg { .. } => 7, + Compression::LossyJpeg { .. } => 34892, + }]), + ), (tags::SAMPLES_PER_PIXEL, Value::Short(vec![o.samples_per_pixel as u16])), (tags::PLANAR_CONFIGURATION, Value::Short(vec![1])), ]; diff --git a/crates/dng/tests/it/equivalence.rs b/crates/dng/tests/it/equivalence.rs new file mode 100644 index 0000000000..dfcfb44f6f --- /dev/null +++ b/crates/dng/tests/it/equivalence.rs @@ -0,0 +1,192 @@ +//! The faster decoder and pipeline must give exactly what the reference implementation of PR +//! #1068 (`filmcraft_dng::reference`) gives: the same raw samples, and the same f32 bit patterns +//! after development, over every Bayer phase, both demosaic qualities, crops with odd origins and +//! at the borders, active-area offsets, linearisation tables, black-level deltas, lossless-JPEG +//! codes above the stated bit depth, LinearRaw, other CFA patterns and tiny pictures. + +use filmcraft_dng::tiff::tags; +use filmcraft_dng::writer::{Compression, Layout, Options, Value, write}; +use filmcraft_dng::{DemosaicQuality, Dng, ProcessOptions}; + +fn noise(n: usize, max: u32, seed: u32) -> Vec { + let mut s = seed.wrapping_mul(2_654_435_761).wrapping_add(1); + (0..n) + .map(|i| { + s ^= s << 13; + s ^= s >> 17; + s ^= s << 5; + // a ramp with noise, some samples clipped high and low + let v = (i as u32 * 37) % (max + 1); + (match s % 23 { + 0 => 0, + 1 => max, + _ => (v ^ (s & 63)).min(max), + }) as u16 + }) + .collect() +} + +fn same_bits(a: &[f32], b: &[f32]) -> bool { + a.len() == b.len() && a.iter().zip(b).all(|(x, y)| x.to_bits() == y.to_bits()) +} + +fn colour_tags() -> Vec<(u16, Value)> { + let cm = [0.9f64, -0.3, -0.05, -0.4, 1.2, 0.2, -0.05, 0.15, 0.6]; + vec![ + (tags::COLOR_MATRIX_1, Value::srationals(&cm)), + (tags::CALIBRATION_ILLUMINANT_1, Value::Short(vec![21])), + (tags::AS_SHOT_NEUTRAL, Value::rationals(&[0.55, 1.0, 0.7])), + ] +} + +/// Raw samples and every development option of `f` agree with the reference. +fn check(f: &[u8], what: &str) { + let d = Dng::parse(f).unwrap(); + assert_eq!(d.raw().unwrap(), d.reference_raw().unwrap(), "{what}: raw"); + for demosaic in [DemosaicQuality::High, DemosaicQuality::Fast] { + let o = ProcessOptions { demosaic, ..Default::default() }; + let (a, b) = (d.develop_with(&o).unwrap(), d.reference_develop(&o).unwrap()); + assert_eq!((a.width, a.height), (b.width, b.height), "{what} {demosaic:?}"); + assert!(same_bits(&a.rgba, &b.rgba), "{what} {demosaic:?}: developed values differ"); + } +} + +#[test] +fn bayer_phases_crops_active_areas_linearisation_and_deltas() { + let phases = [[0u8, 1, 1, 2], [1, 0, 2, 1], [1, 2, 0, 1], [2, 1, 1, 0]]; + let sizes = [(1, 1), (2, 2), (3, 3), (4, 5), (5, 3), (7, 6), (26, 18)]; + for (pi, phase) in phases.iter().enumerate() { + for &(w, h) in &sizes { + for variant in 0..6 { + let px = noise(w * h, 4095, (pi * 100 + w * 7 + h + variant) as u32); + let mut raw_tags = vec![(tags::BLACK_LEVEL, Value::Short(vec![60])), (tags::WHITE_LEVEL, Value::Short(vec![4000]))]; + match variant { + // an active area with an odd offset, a crop with an odd origin inside it + 1 if w > 4 && h > 4 => { + raw_tags.push((tags::ACTIVE_AREA, Value::Long(vec![1, 1, h as u32 - 1, w as u32]))); + raw_tags.push((tags::DEFAULT_CROP_ORIGIN, Value::Long(vec![1, 1]))); + raw_tags.push((tags::DEFAULT_CROP_SIZE, Value::Long(vec![w as u32 - 3, h as u32 - 3]))); + } + // a crop touching the right and bottom borders + 2 if w > 2 && h > 2 => { + raw_tags.push((tags::DEFAULT_CROP_ORIGIN, Value::Long(vec![1, 2.min(h as u32 - 1)]))); + raw_tags.push((tags::DEFAULT_CROP_SIZE, Value::Long(vec![w as u32, h as u32]))); + } + // a linearisation table + 3 => raw_tags.push((tags::LINEARIZATION_TABLE, Value::Short((0..4096u32).map(|v| (v * v / 4096) as u16).collect()))), + // per-site black levels with row and column deltas (the direct path) + 4 => { + raw_tags.push((tags::BLACK_LEVEL_REPEAT_DIM, Value::Short(vec![2, 2]))); + raw_tags.push((tags::BLACK_LEVEL, Value::Short(vec![50, 61, 72, 83]))); + raw_tags.push((tags::BLACK_LEVEL_DELTA_H, Value::srationals(&(0..w).map(|x| x as f64 * 0.5).collect::>()))); + raw_tags.push((tags::BLACK_LEVEL_DELTA_V, Value::srationals(&(0..h).map(|y| -(y as f64) * 0.25).collect::>()))); + } + // per-site black levels without deltas (the lookup table path) + 5 => { + raw_tags.push((tags::BLACK_LEVEL_REPEAT_DIM, Value::Short(vec![2, 2]))); + raw_tags.push((tags::BLACK_LEVEL, Value::Short(vec![50, 61, 72, 83]))); + } + _ => {} + } + let o = Options { + width: w, + height: h, + bits: 12, + cfa: Some((2, 2, phase.to_vec())), + compression: Compression::LosslessJpeg { predictor: 1 + (variant % 7) as u8, components: 1, restart_lines: 0, pt: 0 }, + main_tags: colour_tags(), + raw_tags, + ..Default::default() + }; + check(&write(&px, &o), &format!("phase {phase:?} {w}x{h} variant {variant}")); + } + } + } +} + +#[test] +fn lossless_jpeg_codes_above_the_bit_depth_and_layouts() { + // 16-bit lossless JPEG data in a file that says 12 bits: codes above 4095 reach the decoder + let (w, h) = (40, 22); + let px = noise(w * h, 65_535, 9); + for (layout, components, restart_lines, predictor) in [ + (Layout::Strips { rows: h }, 2, 0, 1), + (Layout::Strips { rows: 5 }, 1, 2, 6), + (Layout::Tiles { width: 16, height: 16 }, 2, 0, 4), + (Layout::Tiles { width: 32, height: 8 }, 4, 1, 7), + ] { + for be in [false, true] { + let o = Options { + width: w, + height: h, + bits: 16, + big_endian: be, + layout, + compression: Compression::LosslessJpeg { predictor, components, restart_lines, pt: 0 }, + main_tags: colour_tags(), + raw_tags: vec![(tags::BITS_PER_SAMPLE, Value::Short(vec![12]))], + ..Default::default() + }; + check(&write(&px, &o), &format!("{layout:?} {components} components restart {restart_lines} be {be}")); + } + } +} + +#[test] +fn uncompressed_linear_raw_mono_and_other_patterns() { + let (w, h) = (19, 11); + for bits in [8, 10, 14, 16] { + let max = (1u32 << bits) - 1; + // three-sample LinearRaw + let px = noise(w * h * 3, max, bits); + let o = Options { + width: w, + height: h, + samples_per_pixel: 3, + bits, + cfa: None, + main_tags: colour_tags(), + layout: Layout::Strips { rows: 4 }, + ..Default::default() + }; + check(&write(&px, &o), &format!("LinearRaw {bits}-bit")); + // monochrome + let px = noise(w * h, max, bits + 1); + let o = Options { width: w, height: h, bits, cfa: None, layout: Layout::Tiles { width: 16, height: 16 }, ..Default::default() }; + check(&write(&px, &o), &format!("mono {bits}-bit")); + } + // non-Bayer colour filter arrays (bilinear for both qualities) + for (r, c, p) in [(1, 3, vec![0u8, 1, 2]), (2, 4, vec![0, 1, 2, 1, 1, 2, 1, 0]), (2, 2, vec![0, 1, 2, 1])] { + let px = noise(w * h, 4095, 3); + let o = Options { width: w, height: h, bits: 12, cfa: Some((r, c, p.clone())), main_tags: colour_tags(), ..Default::default() }; + check(&write(&px, &o), &format!("CFA {r}x{c} {p:?}")); + } +} + +#[test] +fn reduced_resolution_covers_the_target_and_keeps_flat_colours() { + let (w, h) = (64, 40); + // a flat colour: every 2 × 2 site has its own constant level + let pat = [0usize, 1, 1, 2]; + let level = [1000u16, 2000, 2000, 3000]; + let px: Vec = (0..w * h).map(|i| level[pat[(i / w % 2) * 2 + i % 2]]).collect(); + let o = Options { width: w, height: h, bits: 12, main_tags: colour_tags(), ..Default::default() }; + let f = write(&px, &o); + let d = Dng::parse(&f).unwrap(); + let full = d.develop().unwrap(); + for (target, want) in [((64, 40), (64, 40)), ((32, 20), (32, 20)), ((30, 10), (32, 20)), ((16, 10), (16, 10)), ((1, 1), (8, 5)), ((40, 40), (64, 40))] { + let img = d.develop_with(&ProcessOptions { target_size: Some(target), ..Default::default() }).unwrap(); + assert_eq!((img.width, img.height), want, "target {target:?}"); + assert_eq!(img.rgba.len(), want.0 * want.1 * 4); + // a flat field develops to the same colour at every size + for (a, b) in img.rgba.iter().zip(&full.rgba) { + assert!((a - b).abs() < 1e-5, "target {target:?}: {a} vs {b}"); + } + } + // LinearRaw and monochrome reduce too + let lin: Vec = (0..w * h).flat_map(|_| [500u16, 1500, 2500]).collect(); + let o = Options { width: w, height: h, samples_per_pixel: 3, bits: 12, cfa: None, ..Default::default() }; + let img = filmcraft_dng::decode_with(&write(&lin, &o), &ProcessOptions { target_size: Some((10, 10)), ..Default::default() }).unwrap(); + assert_eq!((img.width, img.height), (16, 10)); + assert!((img.rgba[0] - 500.0 / 4095.0).abs() < 1e-6); +} diff --git a/crates/dng/tests/it/lossy.rs b/crates/dng/tests/it/lossy.rs new file mode 100644 index 0000000000..d0d6edfbb3 --- /dev/null +++ b/crates/dng/tests/it/lossy.rs @@ -0,0 +1,130 @@ +//! Lossy JPEG DNG (Compression 34892, DNG 1.4): 8-bit LinearRaw tiles or strips of one or three +//! components, usually with a LinearizationTable taking the 8-bit codes back to linear values. + +use filmcraft_dng::tiff::tags; +use filmcraft_dng::writer::{Compression, Layout, Options, Value, lossy_jpeg, write}; +use filmcraft_dng::{Dng, Error, decode, probe}; + +/// A smooth 8-bit picture (what JPEG is good at). +fn smooth(w: usize, h: usize, spp: usize) -> Vec { + (0..w * h * spp) + .map(|i| { + let (x, y, s) = ((i / spp) % w, (i / spp) / w, i % spp); + (40.0 + 80.0 * (1.0 + ((x as f32 + 3.0 * s as f32) / 9.0).sin() * ((y as f32) / 7.0).cos())) as u16 + }) + .collect() +} + +fn lossy(px: &[u16], w: usize, h: usize, spp: usize, layout: Layout, quality: u8, raw_tags: Vec<(u16, Value)>) -> Vec { + let o = Options { + width: w, + height: h, + samples_per_pixel: spp, + bits: 8, + cfa: None, + layout, + compression: Compression::LossyJpeg { quality }, + raw_tags, + ..Default::default() + }; + write(px, &o) +} + +#[test] +fn lossy_jpeg_tiles_and_strips_one_and_three_components() { + let (w, h) = (40, 28); + for spp in [1, 3] { + let px = smooth(w, h, spp); + for layout in [Layout::Strips { rows: h }, Layout::Strips { rows: 9 }, Layout::Tiles { width: 16, height: 16 }] { + let f = lossy(&px, w, h, spp, layout, 100, vec![]); + let info = probe(&f).unwrap(); + assert_eq!((info.compression, info.bits, info.cfa, info.samples_per_pixel), ("lossy JPEG", 8, false, spp)); + let raw = Dng::parse(&f).unwrap().raw().unwrap(); + assert_eq!(raw.data.len(), px.len()); + // quality 100 is close to the original (colour goes through YCbCr) + let worst = raw.data.iter().zip(&px).map(|(a, b)| (i32::from(*a) - i32::from(*b)).abs()).max().unwrap(); + assert!(worst <= 4, "spp {spp} {layout:?}: off by {worst}"); + // a picture comes out + let img = decode(&f).unwrap(); + assert_eq!((img.width, img.height), (w, h)); + } + } +} + +#[test] +fn samples_are_the_jpeg_decoders_exactly() { + // one strip = one JPEG: our samples are the decoded JPEG's, in raster order + let (w, h) = (24, 16); + for spp in [1, 3] { + let px = smooth(w, h, spp); + let f = lossy(&px, w, h, spp, Layout::Strips { rows: h }, 75, vec![]); + let jpeg = lossy_jpeg(&px, w, h, spp, 75); + let want: Vec = match image::load_from_memory(&jpeg).unwrap() { + image::DynamicImage::ImageLuma8(g) => g.into_raw().into_iter().map(u16::from).collect(), + image::DynamicImage::ImageRgb8(c) => c.into_raw().into_iter().map(u16::from).collect(), + other => panic!("{:?}", other.color()), + }; + assert_eq!(Dng::parse(&f).unwrap().raw().unwrap().data, want, "spp {spp}"); + } +} + +#[test] +fn linearisation_table_expands_the_codes() { + // grey 8-bit codes through a square-law table to 16 bits, as lossy DNG writers do + let (w, h) = (16, 16); + let px = vec![128u16; w * h]; + let table: Vec = (0..256u32).map(|v| (v * v * 65_535 / (255 * 255)) as u16).collect(); + let f = lossy( + &px, + w, + h, + 1, + Layout::Strips { rows: h }, + 90, + vec![(tags::LINEARIZATION_TABLE, Value::Short(table.clone())), (tags::WHITE_LEVEL, Value::Long(vec![65_535]))], + ); + let img = decode(&f).unwrap(); + let code = Dng::parse(&f).unwrap().raw().unwrap().data[0]; + let want = f32::from(table[usize::from(code)]) / 65_535.0; + assert!((127..=129).contains(&code), "{code}"); + assert!((img.rgba[0] - want).abs() < 1e-6 && (img.rgba[0] - 0.252).abs() < 0.01, "{}", img.rgba[0]); +} + +#[test] +fn only_8_bit_linear_raw_and_hostile_lossy_files() { + let (w, h) = (16, 16); + // CFA or 12-bit lossy JPEG is outside the specification + let o = Options { width: w, height: h, bits: 12, raw_tags: vec![(tags::COMPRESSION, Value::Short(vec![34892]))], ..Default::default() }; + let e = decode(&write(&vec![100; w * h], &o)).err(); + assert!(matches!(e, Some(Error::Unsupported(ref m)) if m.contains("lossy JPEG") && m.contains("CFA")), "{e:?}"); + let o = Options { width: w, height: h, bits: 12, cfa: None, raw_tags: vec![(tags::COMPRESSION, Value::Short(vec![34892]))], ..Default::default() }; + assert!(matches!(decode(&write(&vec![100; w * h], &o)), Err(Error::Unsupported(m)) if m.contains("12-bit"))); + // a tiny file claiming a huge picture is refused before allocating it + let mut big = lossy(&smooth(w, h, 3), w, h, 3, Layout::Strips { rows: h }, 80, vec![]); + let ifd = u32::from_le_bytes(big[4..8].try_into().unwrap()) as usize; + let n = u16::from_le_bytes(big[ifd..ifd + 2].try_into().unwrap()) as usize; + for e in (0..n).map(|i| ifd + 2 + 12 * i) { + let tag = u16::from_le_bytes(big[e..e + 2].try_into().unwrap()); + if tag == tags::IMAGE_WIDTH || tag == tags::IMAGE_LENGTH || tag == tags::ROWS_PER_STRIP { + big[e + 8..e + 12].copy_from_slice(&60_000u32.to_le_bytes()); + } + } + assert!(decode(&big).is_err()); + // truncations and mutations + let f = lossy(&smooth(w, h, 3), w, h, 3, Layout::Tiles { width: 16, height: 8 }, 85, vec![]); + for n in 0..f.len() { + assert!(std::panic::catch_unwind(|| decode(&f[..n]).map(|_| ())).is_ok(), "truncated to {n}"); + } + let mut s = 4321u32; + for i in 0..3000 { + let mut g = f.clone(); + for _ in 0..1 + i % 4 { + s ^= s << 13; + s ^= s >> 17; + s ^= s << 5; + let at = s as usize % g.len(); + g[at] = if i % 3 == 0 { 0xFF } else { g[at] ^ (1 << (s >> 29)) }; + } + assert!(std::panic::catch_unwind(|| decode(&g).map(|_| ())).is_ok(), "mutation {i}"); + } +} diff --git a/crates/dng/tests/it/main.rs b/crates/dng/tests/it/main.rs index fd498d961d..d2475f3be5 100644 --- a/crates/dng/tests/it/main.rs +++ b/crates/dng/tests/it/main.rs @@ -1,2 +1,5 @@ +mod equivalence; +mod lossy; +mod opcodes; mod oracle_libraw; mod roundtrip; diff --git a/crates/dng/tests/it/opcodes.rs b/crates/dng/tests/it/opcodes.rs new file mode 100644 index 0000000000..7e019eb4ce --- /dev/null +++ b/crates/dng/tests/it/opcodes.rs @@ -0,0 +1,198 @@ +//! Opcode lists in whole files: each list runs at its stage (OpcodeList1 on the stored samples, +//! OpcodeList2 on linear values before white balance, OpcodeList3 after the demosaic), unknown +//! opcodes follow their optional flag, preview-only opcodes are skipped on the reduced path, and +//! hostile lists are errors, never panics. + +use filmcraft_dng::opcode::encode; +use filmcraft_dng::tiff::tags; +use filmcraft_dng::writer::{Compression, Options, Value, write}; +use filmcraft_dng::{Dng, Error, ProcessOptions, decode_with}; + +/// Develop with the opcode lists applied (they are off by default). +fn decode(bytes: &[u8]) -> Result { + decode_with(bytes, &ProcessOptions { apply_opcodes: true, ..Default::default() }) +} + +fn u(v: &[u32]) -> Vec { + v.iter().flat_map(|x| x.to_be_bytes()).collect() +} + +fn d(v: &[f64]) -> Vec { + v.iter().flat_map(|x| x.to_be_bytes()).collect() +} + +fn f(v: &[f32]) -> Vec { + v.iter().flat_map(|x| x.to_be_bytes()).collect() +} + +fn picture(w: usize, h: usize, seed: u32) -> Vec { + let mut s = seed.wrapping_mul(2_654_435_761).wrapping_add(1); + (0..w * h) + .map(|i| { + s ^= s << 13; + s ^= s >> 17; + s ^= s << 5; + (500 + (i * 13) % 2500 + (s % 64) as usize) as u16 + }) + .collect() +} + +/// A 12-bit RGGB file, white balanced by (0.5, 1, 0.8), with `lists` as (tag, list) pairs. +fn bayer(px: &[u16], w: usize, h: usize, lists: &[(u16, Vec)]) -> Vec { + let mut raw_tags = vec![(tags::BLACK_LEVEL, Value::Short(vec![100])), (tags::WHITE_LEVEL, Value::Short(vec![4095]))]; + raw_tags.extend(lists.iter().map(|(t, l)| (*t, Value::Undefined(l.clone())))); + let o = Options { + width: w, + height: h, + bits: 12, + compression: Compression::LosslessJpeg { predictor: 1, components: 1, restart_lines: 0, pt: 0 }, + main_tags: vec![(tags::AS_SHOT_NEUTRAL, Value::rationals(&[0.5, 1.0, 0.8]))], + raw_tags, + ..Default::default() + }; + write(px, &o) +} + +fn same_bits(a: &[f32], b: &[f32]) -> bool { + a.len() == b.len() && a.iter().zip(b).all(|(x, y)| x.to_bits() == y.to_bits()) +} + +#[test] +fn opcode_list_1_runs_on_the_stored_samples() { + let (w, h) = (16, 12); + let px = picture(w, h, 1); + // MapTable halving every code over the whole raw image = a file holding the halved codes + let table: Vec = (0..4096u32).flat_map(|v| ((v / 2) as u16).to_be_bytes()).collect(); + let list = encode(&[(7, 0, [u(&[0, 0, h as u32, w as u32, 0, 1, 1, 1, 4096]), table].concat())]); + let halved: Vec = px.iter().map(|v| v / 2).collect(); + let a = decode(&bayer(&px, w, h, &[(tags::OPCODE_LIST_1, list)])).unwrap(); + let b = decode(&bayer(&halved, w, h, &[])).unwrap(); + assert!(same_bits(&a.rgba, &b.rgba)); + // raw() stays as stored + let f = bayer(&px, w, h, &[(tags::OPCODE_LIST_1, encode(&[(7, 0, vec![])]))]); + assert_eq!(Dng::parse(&f).unwrap().raw().unwrap().data, px); + assert!(decode(&f).is_err(), "a mandatory MapTable without parameters is malformed"); +} + +#[test] +fn opcode_list_2_runs_before_white_balance() { + let (w, h) = (16, 12); + let px = picture(w, h, 2); + // ScalePerRow 0.5 on every row of the linear values = halving (code − black) + let list = encode(&[(12, 0, [u(&[0, 0, h as u32, w as u32, 0, 1, 1, 1, h as u32]), f(&vec![0.5; h])].concat())]); + let a = decode(&bayer(&px, w, h, &[(tags::OPCODE_LIST_2, list)])).unwrap(); + let halved: Vec = px.iter().map(|&v| ((f64::from(v) - 100.0) / 2.0 + 100.0) as u16).collect(); + let b = decode(&bayer(&halved, w, h, &[])).unwrap(); + for (x, y) in a.rgba.iter().zip(&b.rgba) { + assert!((x - y).abs() < 2e-3, "{x} {y}"); + } + // a gain only on the green sites of the mosaic (row / column pitch 2 from (0, 1) and (1, 0)) + let green = encode(&[ + (12, 0, [u(&[0, 1, h as u32, w as u32, 0, 1, 2, 2, (h / 2) as u32]), f(&vec![0.0; h / 2])].concat()), + (12, 0, [u(&[1, 0, h as u32, w as u32, 0, 1, 2, 2, (h / 2) as u32]), f(&vec![0.0; h / 2])].concat()), + ]); + let g = decode(&bayer(&px, w, h, &[(tags::OPCODE_LIST_2, green)])).unwrap(); + let plain = decode(&bayer(&px, w, h, &[])).unwrap(); + let mean_green = |img: &[f32]| img.chunks(4).map(|p| p[1]).sum::(); + assert!(mean_green(&g.rgba) < 0.1 * mean_green(&plain.rgba), "green sites zeroed"); + // red and blue sites untouched: the red of a red site is the same + assert_eq!(g.rgba[0], plain.rgba[0]); +} + +#[test] +fn opcode_list_3_runs_after_the_demosaic_and_may_trim() { + let (w, h) = (16, 12); + let px = picture(w, h, 3); + let plain = decode(&bayer(&px, w, h, &[])).unwrap(); + // TrimBounds to (2, 1)–(10, 9) on the demosaiced picture + let trim = encode(&[(6, 0, u(&[1, 2, 9, 10]))]); + let f = bayer(&px, w, h, &[(tags::OPCODE_LIST_3, trim)]); + let t = decode(&f).unwrap(); + assert_eq!((t.width, t.height), (8, 8)); + assert_eq!(Dng::parse(&f).unwrap().info().unwrap().width, w, "the info keeps the stored crop"); + for y in 0..8 { + for x in 0..8 { + let a = &t.rgba[(y * 8 + x) * 4..][..4]; + let b = &plain.rgba[((y + 1) * w + x + 2) * 4..][..4]; + for k in 0..3 { + assert!((a[k] - b[k]).abs() < 1e-5, "({x},{y}) {k}: {} {}", a[k], b[k]); + } + } + } + // a warp on the demosaiced picture is allowed (identity changes nothing beyond clipping) + let warp = encode(&[(1, 0, [u(&[1]), d(&[1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5])].concat())]); + let wpd = decode(&bayer(&px, w, h, &[(tags::OPCODE_LIST_3, warp.clone())])).unwrap(); + for (a, b) in wpd.rgba.iter().zip(&plain.rgba) { + assert!((a - b.min(1.0 / 0.5)).abs() < 1e-4, "{a} {b}"); + } + // the same warp before the demosaic is refused (mandatory) or skipped (optional) + let e = decode(&bayer(&px, w, h, &[(tags::OPCODE_LIST_2, warp)])); + assert!(matches!(e, Err(Error::Unsupported(_))), "{:?}", e.err()); + let optional = encode(&[(1, 1, [u(&[1]), d(&[0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5])].concat())]); + let o = decode(&bayer(&px, w, h, &[(tags::OPCODE_LIST_2, optional)])).unwrap(); + assert!(same_bits(&o.rgba, &plain.rgba)); +} + +#[test] +fn unknown_opcodes_and_the_preview_path() { + let (w, h) = (32, 24); + let px = picture(w, h, 4); + let plain = decode(&bayer(&px, w, h, &[])).unwrap(); + // an unknown optional opcode is skipped, a mandatory one names itself + let f = bayer(&px, w, h, &[(tags::OPCODE_LIST_1, encode(&[(4242, 1, vec![1, 2])]))]); + assert!(same_bits(&decode(&f).unwrap().rgba, &plain.rgba)); + let f = bayer(&px, w, h, &[(tags::OPCODE_LIST_3, encode(&[(4242, 0, vec![])]))]); + assert!(matches!(decode(&f), Err(Error::Unsupported(m)) if m.contains("4242"))); + // a preview-only vignette in OpcodeList3 is skipped at half size and applied at full size + let vig = |flags| encode(&[(3, flags, d(&[0.5, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5]))]); + let f = bayer(&px, w, h, &[(tags::OPCODE_LIST_3, vig(2))]); + let half = decode_with(&f, &ProcessOptions { target_size: Some((16, 12)), apply_opcodes: true, ..Default::default() }).unwrap(); + assert_eq!((half.width, half.height), (16, 12)); + let half_plain = decode_with(&bayer(&px, w, h, &[]), &ProcessOptions { target_size: Some((16, 12)), apply_opcodes: true, ..Default::default() }).unwrap(); + assert!(same_bits(&half.rgba, &half_plain.rgba)); + let full = decode(&f).unwrap(); + assert!(full.rgba[0] > plain.rgba[0], "the corner brightened"); + // a mandatory OpcodeList3 opcode keeps full resolution + let f = bayer(&px, w, h, &[(tags::OPCODE_LIST_3, vig(0))]); + let img = decode_with(&f, &ProcessOptions { target_size: Some((16, 12)), apply_opcodes: true, ..Default::default() }).unwrap(); + assert_eq!((img.width, img.height), (w, h)); +} + +#[test] +fn hostile_opcode_lists_are_errors_not_panics() { + let (w, h) = (16, 12); + let px = picture(w, h, 5); + let lists = [ + (tags::OPCODE_LIST_1, encode(&[(5, 0, u(&[0, 2, 1, 3, 3, 0, 0, 4, 4, 6, 6])), (4, 1, u(&[0, 0]))])), + (tags::OPCODE_LIST_2, encode(&[(9, 0, [u(&[0, 0, 12, 16, 0, 1, 1, 1, 2, 2]), d(&[1.0, 1.0, 0.0, 0.0]), u(&[1]), f(&[1.0, 1.1, 1.2, 1.3])].concat())])), + (tags::OPCODE_LIST_3, encode(&[(3, 0, d(&[0.1, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5])), (6, 0, u(&[1, 1, 11, 15]))])), + ]; + let f = bayer(&px, w, h, &lists); + if let Err(e) = decode(&f) { + panic!("{e}"); + } + let mut s = 777u32; + for i in 0..3000 { + let mut g = f.clone(); + for _ in 0..1 + i % 3 { + s ^= s << 13; + s ^= s >> 17; + s ^= s << 5; + let at = s as usize % g.len(); + g[at] = if i % 4 == 0 { 0xFF } else { g[at] ^ (1 << (s >> 29)) }; + } + assert!(std::panic::catch_unwind(|| decode(&g).map(|_| ())).is_ok(), "mutation {i}"); + } +} + +#[test] +fn opcode_lists_are_ignored_unless_asked_for() { + let (w, h) = (16, 12); + let px = picture(w, h, 7); + let table: Vec = (0..4096u32).flat_map(|v| ((v / 2) as u16).to_be_bytes()).collect(); + let list = encode(&[(7, 0, [u(&[0, 0, h as u32, w as u32, 0, 1, 1, 1, 4096]), table].concat())]); + let f = bayer(&px, w, h, &[(tags::OPCODE_LIST_1, list)]); + let plain = filmcraft_dng::decode(&bayer(&px, w, h, &[])).unwrap(); + assert!(same_bits(&filmcraft_dng::decode(&f).unwrap().rgba, &plain.rgba), "default ignores the lists"); + assert!(!same_bits(&decode(&f).unwrap().rgba, &plain.rgba), "apply_opcodes uses them"); +} diff --git a/crates/dng/tests/it/roundtrip.rs b/crates/dng/tests/it/roundtrip.rs index af9d8e287f..d586c04ed9 100644 --- a/crates/dng/tests/it/roundtrip.rs +++ b/crates/dng/tests/it/roundtrip.rs @@ -5,7 +5,7 @@ use filmcraft_dng::color::{D50, D65, bradford, inverse, mul, rec709_to_xyz, xy_to_xyz}; use filmcraft_dng::tiff::tags; use filmcraft_dng::writer::{Compression, Layout, Lj92Params, Options, Value, lj92, write}; -use filmcraft_dng::{Dng, Error, Timecode, decode, ljpeg, probe, sniff}; +use filmcraft_dng::{DemosaicQuality, Dng, Error, Highlights, ProcessOptions, Timecode, WhiteBalance, decode, decode_with, ljpeg, probe, sniff}; /// A deterministic pseudo-random sequence. fn noise(n: usize, bits: u32, seed: u32) -> Vec { @@ -365,7 +365,88 @@ fn cyclic_ifds_huge_sizes_and_bad_offsets() { assert!(matches!(decode(&write(&px, &o)), Err(Error::Malformed(_)))); // unsupported compressions and depths are named let o = Options { raw_tags: vec![(tags::COMPRESSION, Value::Short(vec![34892]))], ..Default::default() }; - assert_eq!(decode(&write(&px, &o)).err(), Some(Error::Unsupported("lossy JPEG DNG".into()))); + assert!(matches!(decode(&write(&px, &o)), Err(Error::Unsupported(m)) if m.contains("lossy JPEG"))); +} + +#[test] +fn process_options_default_is_the_files_development() { + let (w, h) = (24, 16); + let px = picture(w, h, 1, 12, 11); + // a camera that is not Rec. 709 already, white balanced away from D65 + let cm = [[0.9, -0.3, -0.05], [-0.4, 1.2, 0.2], [-0.05, 0.15, 0.6]]; + let main = vec![ + (tags::COLOR_MATRIX_1, Value::srationals(&cm.iter().flatten().copied().collect::>())), + (tags::CALIBRATION_ILLUMINANT_1, Value::Short(vec![21])), + (tags::AS_SHOT_NEUTRAL, Value::rationals(&[0.55, 1.0, 0.7])), + ]; + let o = Options { width: w, height: h, bits: 12, main_tags: main, raw_tags: vec![(tags::WHITE_LEVEL, Value::Short(vec![4095]))], ..Default::default() }; + let f = write(&px, &o); + let base = decode(&f).unwrap(); + let same = decode_with(&f, &ProcessOptions::default()).unwrap(); + assert_eq!(base.rgba, same.rgba, "default options change nothing"); + // highlight recovery is not implemented yet: identical to clipping + let rec = decode_with(&f, &ProcessOptions { highlights: Highlights::Recover, ..Default::default() }).unwrap(); + assert_eq!(base.rgba, rec.rgba); + // +1 stop doubles every value + let up = decode_with(&f, &ProcessOptions { exposure: 1.0, ..Default::default() }).unwrap(); + for (a, b) in base.rgba.iter().zip(&up.rgba).filter(|(a, _)| **a != 1.0) { + assert!((b - 2.0 * a).abs() <= 1e-5 * a.max(1.0), "{a} {b}"); + } + let nan = decode_with(&f, &ProcessOptions { exposure: f64::NAN, ..Default::default() }).unwrap(); + assert_eq!(base.rgba, nan.rgba, "a non-finite exposure is ignored"); + // the as-shot white as temperature / tint develops like the as-shot neutral + let (t, tint) = Dng::parse(&f).unwrap().as_shot_temperature_tint().unwrap().unwrap(); + assert!((1_000.0..15_000.0).contains(&t), "{t}"); + let custom = decode_with(&f, &ProcessOptions { white_balance: WhiteBalance::Custom { temperature: t, tint }, ..Default::default() }).unwrap(); + for (a, b) in base.rgba.iter().zip(&custom.rgba) { + assert!((a - b).abs() < 2e-3, "{a} {b}"); + } + // a warmer setting makes the picture bluer (the white balance compensates a warmer light) + let warm = decode_with(&f, &ProcessOptions { white_balance: WhiteBalance::Custom { temperature: t * 0.6, tint }, ..Default::default() }).unwrap(); + let ratio = |img: &[f32]| img.chunks(4).map(|p| p[2] / p[0].max(1e-6)).sum::(); + assert!(ratio(&warm.rgba) > ratio(&base.rgba)); + assert!(decode_with(&f, &ProcessOptions { white_balance: WhiteBalance::Custom { temperature: f64::NAN, tint: 0.0 }, ..Default::default() }).is_err()); + // the fast demosaic keeps flat fields flat and differs from the high-quality one on detail + let fast = decode_with(&f, &ProcessOptions { demosaic: DemosaicQuality::Fast, ..Default::default() }).unwrap(); + assert_eq!((fast.width, fast.height), (w, h)); + assert_ne!(fast.rgba, base.rgba); +} + +#[test] +fn process_options_serde_round_trip() { + let o = ProcessOptions { + white_balance: WhiteBalance::Custom { temperature: 5600.0, tint: -12.5 }, + exposure: 0.75, + highlights: Highlights::Recover, + demosaic: DemosaicQuality::Fast, + target_size: Some((1920, 1080)), + ..Default::default() + }; + let j = serde_json::to_string(&o).unwrap(); + assert_eq!(serde_json::from_str::(&j).unwrap(), o); + // missing fields take their defaults + assert_eq!(serde_json::from_str::("{}").unwrap(), ProcessOptions::default()); +} + +#[test] +fn positive_tint_makes_a_grey_more_magenta() { + // a grey under D65 seen by a camera whose white-balanced values are Rec. 709 + let (w, h) = (16, 12); + let px = vec![2000u16; w * h]; + let mut main = rec709_camera_tags([1.0, 1.0, 1.0]); + main.retain(|t| t.0 != tags::FORWARD_MATRIX_1); + let o = Options { width: w, height: h, bits: 12, main_tags: main, raw_tags: vec![(tags::WHITE_LEVEL, Value::Short(vec![4095]))], ..Default::default() }; + let f = write(&px, &o); + let (t, tint) = Dng::parse(&f).unwrap().as_shot_temperature_tint().unwrap().unwrap(); + assert!((t - 6504.0).abs() < 60.0 && (8.0..12.0).contains(&tint), "as shot D65: {t} K tint {tint}"); + let dev = |tint: f64| { + let img = decode_with(&f, &ProcessOptions { white_balance: WhiteBalance::Custom { temperature: t, tint }, ..Default::default() }).unwrap(); + let p = &img.rgba[(6 * w + 6) * 4..][..3]; + p[1] / (0.5 * (p[0] + p[2])) + }; + let (minus, zero, plus) = (dev(tint - 50.0), dev(tint), dev(tint + 50.0)); + assert!((zero - 1.0).abs() < 0.01, "as-shot grey stays grey: {zero}"); + assert!(plus < zero && zero < minus, "green over red/blue: {minus} {zero} {plus}"); } #[test]