From 0da3928eee4880139488e28ec6f444ee2dd722fc Mon Sep 17 00:00:00 2001 From: "Dylan M. Taylor" Date: Wed, 11 Nov 2026 20:22:37 -0500 Subject: [PATCH] craft-raw: Panasonic RW2, RWL and RAW (#1000) * craft-raw: the file-format layer of camera raw files, with DNG and CRW A new crate for #587 and #778: bytes in, the sensor data out as a neutral RawImage (samples, CFA, levels, active area, crop, orientation, colour tags, opcode lists, profile tables), and back for DNG. One Cargo feature per format: dng (read and write), jxl (JPEG XL DNG), crw (through craft-crw). The base is LightCraft's crates/raw at 5ff567a4: dng, dngwrite, tiffraw, unpack, ljpeg, jxl, semantic, preview, the parsing half of opcodes, profile and gaintable, the model and probe. Demosaicing and everything after it stay with the develop pipeline. Ported from PhotoCraft's and FilmCraft's DNG readers, each with a test: the default crop's per-tag defaults, an unusable crop falling back to the active area, non-finite white levels, the white default of a linearized image, the AsShotWhiteXY range, 8-bit-only lossy JPEG, lossy tiles that must cover their part of the image, and the CinemaDNG frame tags. * craft-raw: Hasselblad 3FR and FFF (feature hasselblad) RawFormat::CfaTiff decodes: the DNG-layout reader plus what these files add. Uncompressed data or Hasselblad's lossless JPEG (selection value 8), RGGB when no CFAPattern is stated, an FFF's colour matrix from the 3FR header it embeds, and the 8-bit RGB preview handed on as a baseline TIFF. Moved from LightCraft's vendor/hasselblad.rs at 5ff567a4 with its tests. Ported from PhotoCraft's reader: four BlackLevel values without a BlackLevelRepeatDim are one per CFA site in a file that is not a DNG. * craft-raw: Minolta MRW, Epson ERF, Mamiya MEF and Kodak DCS520C Four formats only LightCraft decodes, one feature each (mrw, erf, mef, kodak), moved from its vendor modules at 5ff567a4 with their tests; the MRW block container comes from lightcraft-meta. src/shared.rs holds the level helpers the maker decoders share. ERF and MEF state no white balance: the estimate LightCraft made from the preview is develop-side work, so wb_multipliers is None for them here. * craft-raw: shared scan budget for tagged previews, tile-sized lossy JPEG budget, bound on Deflate tile padding Ported from storytold/lightcraft#3695. Candidate preview ranges are deduplicated and charged to one byte budget and one header-walk budget; a lossy JPEG decoded for a tile may hold no more samples than the tile rounded up to whole MCUs (and never more than Limits::max_lossy_tile_samples); a Deflate chunk whose nominal size is mostly padding outside the image is refused before it is inflated. * craft-raw: lossless JPEG sample bound per referenced byte; previews larger than the file are passed over Ported from storytold/lightcraft#3910 (the parts whose code is in craft-raw today). check_image refuses a lossless JPEG image that declares more than 32 samples per referenced byte, counted once per chunk that points at the bytes so identical tiles stored once still decode. The Hasselblad RGB preview search passes over an image whose pixel bytes exceed the file and still takes a smaller complete one. * craft-raw: test for the Deflate padding bound * craft-raw: Pentax PEF and Samsung SRW Moved from LightCraft crates/raw (vendor/pef.rs, vendor/srw.rs) at 5ff567a4. The Huffman difference decoder and bit reader that PEF shares with the Nikon, Panasonic and Sony codings become the huffman module. * craft-raw: Canon CR2, with sRAW and mRAW Moved from LightCraft crates/raw (vendor/cr2.rs, vendor/cr2/sraw.rs) at 5ff567a4. The white-balance presets a CR2 lists come back as data in ColorData::kelvin_presets (the ordering from LightCraft's kelvin_presets.rs; the temperature and tint arithmetic stays with the developer). * craft-raw: an sRAW preview larger than the file is refused before its buffer is reserved Ported from storytold/lightcraft#3910 (cr2::sraw::read_preview). * WIP craft-raw: Nikon NEF/NRW (does not compile yet: limits threading, develop() calls in tests) * craft-raw: Nikon NEF/NRW, with PhotoCraft's maker-note compression routing and 12-bit D4 white level * craft-raw: Panasonic RW2 / RWL / RAW, every raw format * craft-raw: RW2 format 8 refuses strips stating fewer bits than samples Ported from storytold/lightcraft#3910. * craft-raw: RW2 format 8 refuses samples above the declared bit depth Ported from storytold/photocraft#11838. * craft-raw: the RW2 format 8 robustness seed carries a valid stream Its random payload decoded only while out-of-depth samples were clamped; mutations now start from a file that decodes. --------- Co-authored-by: radialmonster Co-authored-by: wnanbei --- NOTICE | 2 +- README.md | 1 + crates/raw/Cargo.toml | 6 +- crates/raw/README.md | 8 + crates/raw/src/lib.rs | 6 + crates/raw/src/rw2.rs | 1348 ++++++++++++++++++++++++++ crates/raw/src/shared.rs | 164 ++++ crates/raw/tests/it/vendor_robust.rs | 114 +++ 8 files changed, 1646 insertions(+), 3 deletions(-) create mode 100644 crates/raw/src/rw2.rs diff --git a/NOTICE b/NOTICE index eb22134a..7421e32d 100644 --- a/NOTICE +++ b/NOTICE @@ -93,7 +93,7 @@ crates/raw (craft-raw) is moved from LightCraft, with behaviour ported from two - LightCraft, crates/raw (src/lib.rs, dng.rs, dngwrite.rs, tiffraw.rs, unpack.rs, ljpeg.rs, jxl.rs, opcodes.rs, profile.rs, gaintable.rs, semantic.rs, preview.rs, vendor/crw.rs, vendor/hasselblad.rs, vendor.rs (shared level helpers), vendor/mrw.rs, vendor/erf.rs, vendor/mef.rs, vendor/kodak.rs, vendor/pef.rs, - vendor/srw.rs, vendor/cr2.rs, vendor/cr2/sraw.rs, vendor/nef.rs, vendor/nefc.rs, vendor/nef_lens.rs, docs/nikon-lens-corrections.md, kelvin_presets.rs (the preset ordering), tests/it/vendor_robust.rs and their tests) and crates/meta/src/mrw.rs (the MRW block container), + vendor/srw.rs, vendor/cr2.rs, vendor/cr2/sraw.rs, vendor/nef.rs, vendor/nefc.rs, vendor/nef_lens.rs, docs/nikon-lens-corrections.md, vendor/rw2.rs, kelvin_presets.rs (the preset ordering), tests/it/vendor_robust.rs and their tests) and crates/meta/src/mrw.rs (the MRW block container), revision 5ff567a4a2e067bef1b6bce152c11e7f727ccbaf. https://github.com/storytold/lightcraft Copyright (c) 2026 ArtCraft Team and the LightCraft contributors. diff --git a/README.md b/README.md index 4292ead6..3a1f14b8 100644 --- a/README.md +++ b/README.md @@ -146,6 +146,7 @@ craft- = { git = "https://github.com/storytold/craft-libs", rev = " Result { RawFormat::Cr2 => cr2::decode(bytes, mode, limits), #[cfg(feature = "nef")] RawFormat::Nef | RawFormat::Nrw => nef::decode(bytes, mode, limits), + #[cfg(feature = "rw2")] + RawFormat::Rw2 => rw2::decode(bytes, mode, limits), #[cfg(feature = "pef")] RawFormat::Pef => pef::decode(bytes, mode, limits), #[cfg(feature = "srw")] diff --git a/crates/raw/src/rw2.rs b/crates/raw/src/rw2.rs new file mode 100644 index 00000000..9cee4481 --- /dev/null +++ b/crates/raw/src/rw2.rs @@ -0,0 +1,1348 @@ +//! Panasonic RW2, Leica RWL and the older Panasonic / Leica `.RAW` files: every raw encoding these cameras write. +//! +//! **Clean-room.** No raw-decoder source code (dcraw, LibRaw, rawspeed, rawloader/rawler, darktable, RawTherapee, +//! libopenraw, ExifTool's Perl code, …) was read or consulted. What we used: +//! +//! - Tag documentation: ExifTool's PanasonicRaw tag-name tables (IFD0 `0x0002`–`0x0007` sensor size and borders, +//! `0x0009` CFAPattern with its four values, `0x000a` BitsPerSample, `0x000b` Compression, `0x000e`–`0x0010` +//! LinearityLimit, `0x0017`/`0x0037` ISO, `0x001c`–`0x001e` black levels, `0x0024`–`0x0026` WB levels, `0x002d` +//! RawFormat, `0x002e` JpgFromRaw, `0x002f`–`0x0032` crop, `0x0118` RawDataOffset) and the libopenraw format +//! notes (prose: the `IIU` magic; Leica RWL files are RW2). +//! - Our own black-box analysis of 178 CC0 files from raw.pixls.us, about 120 bodies from the DMC-LX1 (2005) to +//! the DC-S1RM2 (2025), every encoding below. The measurement that settled each rule is noted in brackets. +//! +//! Encodings, by raw format (tag `0x002d`; the oldest files have none and are told apart by compression `0x000b`): +//! +//! - **16-bit words** (compression 34828 and 34830: DMC-LX1, DMC-FZ50, Digilux 3): little-endian words at +//! `StripOffsets`, the sample left-justified [every value is a multiple of 2^(16 − bits)]. +//! - **Packed** (formats 2, 5 and 7): 16-byte blocks read as a little-endian 128-bit number holding ⌊128 / bits⌋ +//! samples from bit 0 up (12-bit: 10, 14-bit: 9). Format 5 stores every whole 0x4000-byte chunk rotated (its +//! logical content is bytes `0x1ff8..0x4000` followed by `0..0x1ff8`); formats 2 (DMC-FZ8) and 7 (DC-S5 with +//! later firmware) are not rotated [vertical same-colour differences 20–180× smaller with the right choice]. +//! - **Format 4** (most bodies, 2008–2022): chunks rotated as format 5; each 16-byte block, read as a little-endian +//! 128-bit number from bit 127 down, codes 14 pixels of a row. Before pixels 2, 5, 8 and 11 comes a 2-bit scale +//! `s` for the next three pixels, step `m` = 1, 2, 4 or 16 [per-bit statistics show the field boundaries]. The +//! first pixel of each colour (parity of the pixel index `i`) is a seed: 8 bits `h`, then, if `h ≠ 0` or +//! `i ≥ 12`, 4 bits `l`, value `16h + l`. `h = 0` (with `i < 12`) marks a mapped-out defective pixel, decoded as +//! 0, and the next pixel of that colour is read as the seed instead [these blocks decode smoothly only this way, +//! 10× closer to the rows above and below, and they sit at the same sensor sites in every shot of a body]. Every +//! later pixel is an 8-bit code `j` against the colour's previous value `p`: `j = 0` keeps `p` [3–4× closer to +//! the neighbouring rows than reading it as a difference]; otherwise `p + (j − 128)·m` when `s < 3` and +//! `p ≥ 128·m`, else `j·m + (p mod m)` [scale-3 codes track the absolute level, ≈ level / 16; values stay on the +//! camera's 4k + 3 lattice (GX80, 99.9–100 % per branch) only when the absolute reading keeps the low bits of `p`]. +//! - **Format 6** (DC-S1, S1R, S1H, S5; DC-GH5M2): not rotated. 14-bit: two 14-bit seeds, then three groups of a +//! 2-bit scale and three 10-bit codes (11 pixels; the 4 lowest bits are unused); 12-bit: two 12-bit seeds and four +//! groups of a 2-bit scale and three 8-bit codes (14 pixels). Read from bit 127 down and reconstructed as format 4 +//! with the bias 512 (14-bit) or 128 (12-bit) [code statistics per scale as in format 4; flat saturated blocks +//! decode flat only when the absolute reading keeps the low bits]. Seeds are never deferred and zero codes don't +//! occur in any sample; a zero code keeps the previous value as in format 4. +//! - **Format 8** (DC-GH6 and later: GH6, GH7, G9M2, S5M2, S5M2X, S9, S1RM2, S1M2): the image is split into vertical +//! strips listed in `u16` arrays (a count, then the values; 32-bit values as low/high pairs): `0x0044` absolute +//! offsets, `0x0045` left columns, `0x0046` lengths in bits, `0x0047` widths, `0x0048` heights. Each strip is an +//! independent bit stream whose bits are taken from each byte starting at bit 0 [a constant run at a strip's end +//! reads as repeats of one code only in this order; MSB-first reading gives code frequencies of exactly +//! 2^−length, the signature of random bits]. Tag `0x0040` lists (length, code) for 17 prefix codes, the difference +//! categories 0–16 of ITU-T T.81 (category `c`, then `c` additional bits, a leading 0 meaning negative). Pixels +//! are coded two rows at a time, cell by 2×2 cell, column by column (top-left, bottom-left, top-right, +//! bottom-right) [in row order the two greens trade sites: an edge's even and odd rows of green then sit ~2 px +//! apart with red and blue between them]; each is +//! predicted from the same site of the previous cell, the first cell of a row pair from the first cell of the +//! previous pair, 0 at the top [the strips decode to the camera JPEG's image without seams]. The streams end a +//! few padding rows before their stated length. Tags `0x0039`–`0x003f` and `0x0041` carry the same +//! values in every sample; files with other values are reported as unsupported rather than guessed at. The high-resolution +//! mode of the DC-S1RM2 (16296 × 10856, `0x0043` = 2) tiles the image in two rows of four: `0x0045` then lists (x, y) +//! `u16` pairs (y is 0 in single-row files, so the 32-bit reading used for them is the same numbers), the tiles are +//! 5632 and 5224 rows high and are decoded exactly like strips [the tiles' areas add up to the sensor, and `0x0043` equals +//! the number of distinct tile tops in every sample]; the decoder takes the geometry from the lists, not from the body. +//! +//! Levels and geometry: +//! +//! - Black: per colour from `0x001c`–`0x001e`. Formats 4 and 6 store every sample 15 above it [flat dark regions +//! (σ ≈ 1–3) of ~100 files sit 13–20 above the tag, e.g. 143.1 on a D-Lux 7 tagged 128], the other formats +//! don't [131.4 on a format-5 G9 tagged 127.5; 128.7 on a format-8 S9 tagged 128]. Files that record no black +//! (tags 0: bodies up to ~2010) use 15 [their darkest flat regions sit at 15–22]. +//! - White: the saturation plateau of the active area's samples, else the full range of the bit depth. The linearity +//! limits (`0x000e`–`0x0010`) are not used: extended low-ISO shots record limits far below their samples (TZ200D +//! at ISO 100: 3277, samples up to 3529) while the camera's own JPEG renders them against the full range [the +//! JPEG/raw brightness ratio matches normal shots of the same bodies only then], the Digilux 3 records 4095 but +//! saturates at 15871 of 16383, and the GF1 records 16383 for 12-bit samples. +//! - Samples of 0 mark mapped-out defects (format 4 above; isolated zeros between normal neighbours in the 16-bit +//! files) and are interpolated from the same colour two pixels away. +//! - The CFA pattern applies at sample (0, 0) of the data, not at the sensor borders [the two files with odd borders +//! show the greens on the diagonal this predicts]; red and blue were checked against every file's JPEG. +//! - Active area: the sensor borders. Default crop: `0x002f`–`0x0032`, the in-camera aspect ratio (e.g. 4:3 or 1:1 +//! on a 3:2 sensor; the embedded JPEG always shows the whole sensor). Some bodies write zeros there: ignored. +//! +//! Lens distortion (tag `0x0119`, ExifTool name DistortionInfo; ExifTool's names for its words were the only outside +//! input, their meaning below is our own measurement). The tag holds 16 little-endian 16-bit words. The camera's +//! own correction, which its embedded JPEG shows, is applied as an `OpcodeList3` `WarpRectilinear` (the DNG form +//! LightCraft already applies under "Enable Profile Corrections"): +//! +//! - Word 7, low 4 bits: 1 = correction on [0 on the DMC-GH1, whose JPEG matches the uncorrected render]. +//! - Word 12: the normalisation radius in pixels, the half diagonal of the active area [2870 for 4592 × 3448 +//! (2871.2), 3240 for 5184 × 3888, 3611 for 5776 × 4336, 3288 for the FZ1000 II's whole 5472 × 3648 sensor +//! even in its 4:3 mode]. +//! - Words 8, 4 and 11 (ExifTool: DistortionParam08, 04, 11), divided by 32768, are `k2`, `k4`, `k6` of the +//! correction: a sensor point at radius ρ (in units of word 12, about the centre of the active area) appears in +//! the corrected image at ρ·(1 + k2 ρ² + k4 ρ⁴ + k6 ρ⁶). +//! - Word 5 (ExifTool: DistortionScale), divided by 32768, plus 1, is the zoom `z` the camera crops by afterwards: +//! output radius = corrected radius / z [z < 1 when the corrected corners would otherwise be empty]. +//! +//! [These rules were chosen among alternative word orders, powers, directions and readings of the zoom by fitting +//! radial curves between our uncorrected renders and the embedded JPEGs of 38 DMC-GX85 photos (12–32 mm and 25 mm +//! lenses). With no free parameter the words then predict the curves of 35 held-out photos from 12 other bodies and +//! 9 other lenses (9.1–105 mm: two compacts, Panasonic, Leica, Olympus, Sigma and Samyang lenses) within 0.12 % of +//! the half diagonal, against up to 13 % uncorrected; corrected renders of all 79 photos line up with the camera +//! JPEG within 0.75 px at 1440 px at every corner measured.] Words 2, 3, 6, 9, 10 +//! and 13 are not used (word 13 is always word 9 + 500; word 10 is fixed per body); words 0, 1, 14 and 15 change +//! from shot to shot. The correction is radial, so it is fitted by the DNG form's odd polynomial `u·(kr0 + kr1 u² + +//! kr2 u⁴ + kr3 u⁶)` over the output (within 0.4 px on a 20 MP sensor); implausible values (a radius far from the +//! active area's, a zoom beyond ±25 %, a correction that folds over) leave the photo uncorrected. +//! +//! Not supported ([`RawError::Unsupported`], the embedded preview is shown instead): compression 34826, raw formats +//! 1 and 3, and bit depths other than the ones above (no samples of any of them). + +use crate::huffman::{Bits, Huffman, diff}; +use crate::shared::{active_rows, white_from_data}; +use crate::{ + BlackLevel, Camera, Cfa, ColorData, Limits, Mode, Opcode, OpcodeLists, Orientation, RawData, RawError, RawFormat, RawImage, Rect, Result, +}; +use craft_tiff::{ByteOrder, Ifd, Tiff, Value, tags as t}; +use rayon::prelude::*; + +const SENSOR_WIDTH: u16 = 0x0002; +const SENSOR_HEIGHT: u16 = 0x0003; +const BORDER_TOP: u16 = 0x0004; +const BORDER_LEFT: u16 = 0x0005; +const BORDER_BOTTOM: u16 = 0x0006; +const BORDER_RIGHT: u16 = 0x0007; +const CFA_PATTERN: u16 = 0x0009; +const BITS: u16 = 0x000a; +const COMPRESSION: u16 = 0x000b; +const LINEARITY_LIMIT: [u16; 3] = [0x000e, 0x000f, 0x0010]; +const BLACK_RGB: [u16; 3] = [0x001c, 0x001d, 0x001e]; +const WB_RGB: [u16; 3] = [0x0024, 0x0025, 0x0026]; +const RAW_FORMAT: u16 = 0x002d; +/// The in-camera aspect-ratio crop: top, left, bottom, right in sensor coordinates. +const CROP: [u16; 4] = [0x002f, 0x0030, 0x0031, 0x0032]; +const RAW_OFFSET: u16 = 0x0118; +const HUFFMAN_CODES: u16 = 0x0040; +const STRIP_OFFSETS: u16 = 0x0044; +const STRIP_LEFT: u16 = 0x0045; +const STRIP_BITS: u16 = 0x0046; +const STRIP_WIDTHS: u16 = 0x0047; +const STRIP_HEIGHTS: u16 = 0x0048; +/// Format-8 tags of unknown meaning and the values every sample carries (as `u16` arrays). +const FORMAT8_FIXED: [(u16, &[u16]); 8] = [ + (0x0039, &[6, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1]), + (0x003a, &[6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]), + (0x003b, &[65535]), + (0x003c, &[0]), + (0x003d, &[0]), + (0x003e, &[0]), + (0x003f, &[0]), + (0x0041, &[17, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]), +]; +/// Format-8 tile rows (ExifTool: no name): 1 for a single row of strips, 2 for the 16296 × 10856 high-resolution mode. +const TILE_ROWS: u16 = 0x0043; +/// The camera's lens distortion correction (ExifTool name: DistortionInfo); see the module docs. +const DISTORTION_INFO: u16 = 0x0119; +/// Fixed-point one of the distortion coefficients and zoom. +const DISTORTION_ONE: f64 = 32768.0; +/// Largest error of the fitted inverse (relative to the half diagonal) accepted: real lenses fit within 1.2e-4. +const DISTORTION_MAX_FIT_ERROR: f64 = 1e-3; +/// Most strips seen in a file is 8 (the 16296 × 10856 high-resolution mode: 4 columns by 2 rows of tiles). +const MAX_STRIPS: usize = 64; + +const CHUNK: usize = 0x4000; +const ROTATE: usize = 0x1ff8; +/// Formats 4 and 6 store samples this far above the recorded black level. +const CODEC_BLACK_OFFSET: f32 = 15.0; +/// Black level of files that record none. +const DEFAULT_BLACK: f32 = 15.0; + +/// How the samples are stored. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +enum Encoding { + /// Little-endian 16-bit words, samples left-justified. + Words, + /// Continuous 16-byte blocks; `rotated`: whole 0x4000-byte chunks are stored rotated by 0x1ff8. + Blocks { kind: Block, rotated: bool }, + /// Format 8: prefix-coded strips. + Strips, +} + +/// What one 16-byte block holds. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +enum Block { + /// ⌊128 / bits⌋ samples from bit 0 up. + Packed { + bits: u32, + }, + Format4, + Format6 { + bits14: bool, + }, +} + +impl Block { + fn per(self) -> usize { + match self { + Block::Packed { bits } => (128 / bits.max(1)) as usize, + Block::Format4 | Block::Format6 { bits14: false } => 14, + Block::Format6 { bits14: true } => 11, + } + } + + fn decode(self, b: &[u8; 16], out: &mut [u16]) { + match self { + Block::Packed { bits } => packed(b, bits, out), + Block::Format4 => format4(b, out), + Block::Format6 { bits14 } => format6(b, bits14, out), + } + } +} + +fn encoding(ifd: &Ifd, bits: u32) -> Result { + let blocks = |kind, rotated| Ok(Encoding::Blocks { kind, rotated }); + match (ifd.u16(RAW_FORMAT), ifd.u16(COMPRESSION)) { + (Some(2), _) if bits == 12 => blocks(Block::Packed { bits }, false), + (Some(4), _) if bits == 12 => blocks(Block::Format4, true), + (Some(5), _) if bits == 12 || bits == 14 => blocks(Block::Packed { bits }, true), + (Some(6), _) if bits == 12 || bits == 14 => blocks(Block::Format6 { bits14: bits == 14 }, false), + (Some(7), _) if bits == 14 => blocks(Block::Packed { bits }, false), + (Some(8), _) if (12..=16).contains(&bits) => Ok(Encoding::Strips), + (None, Some(34828 | 34830)) if (12..=16).contains(&bits) => Ok(Encoding::Words), + (Some(f), _) => Err(RawError::Unsupported(format!("Panasonic raw format {f} with {bits}-bit samples"))), + (None, c) => Err(RawError::Unsupported(format!("Panasonic raw without a raw format (compression {c:?})"))), + } +} + +/// Bit fields of a 16-byte block taken as a little-endian 128-bit number, read from bit 127 down. +struct Fields { + word: u128, + left: u32, +} + +impl Fields { + fn new(b: &[u8; 16]) -> Fields { + Fields { word: u128::from_le_bytes(*b), left: 128 } + } + + /// The next `n` (1..=16) bits. + #[inline] + fn take(&mut self, n: u32) -> i32 { + self.left = self.left.saturating_sub(n); + ((self.word >> self.left) & ((1u128 << n) - 1)) as i32 + } +} + +/// Step and absolute flag of the 2-bit scales of formats 4 and 6. +const STEPS: [(i32, bool); 4] = [(1, false), (2, false), (4, false), (16, true)]; + +/// Formats 4 and 6: the next value of a colour from its previous value `p` and code `j` (see the module docs). +#[inline] +fn next(p: i32, j: i32, (m, absolute): (i32, bool), bias: i32) -> i32 { + if j == 0 { + p + } else if !absolute && p >= bias * m { + p + (j - bias) * m + } else { + j * m + (p & (m - 1)) + } +} + +#[inline] +fn put(out: &mut [u16], i: usize, v: i32) { + if let Some(o) = out.get_mut(i) { + *o = v.clamp(0, u16::MAX as i32) as u16; + } +} + +fn packed(b: &[u8; 16], bits: u32, out: &mut [u16]) { + let word = u128::from_le_bytes(*b); + let mask = (1u128 << bits.min(16)) - 1; + for (i, o) in out.iter_mut().enumerate().take((128 / bits.max(1)) as usize) { + *o = ((word >> (i as u32 * bits)) & mask) as u16; + } +} + +fn format4(b: &[u8; 16], out: &mut [u16]) { + let mut f = Fields::new(b); + let (mut prev, mut seeded, mut step) = ([0i32; 2], [false; 2], STEPS[0]); + for i in 0..14 { + if i % 3 == 2 { + step = STEPS[f.take(2) as usize & 3]; + } + let c = i & 1; + let v = if seeded[c] { + prev[c] = next(prev[c], f.take(8), step, 128); + prev[c] + } else { + let hi = f.take(8); + if hi != 0 || i >= 12 { + prev[c] = hi << 4 | f.take(4); + seeded[c] = hi != 0; + prev[c] + } else { + 0 // a defect marker; this colour's seed comes with its next pixel + } + }; + put(out, i, v); + } +} + +fn format6(b: &[u8; 16], bits14: bool, out: &mut [u16]) { + let (seed, code, groups, bias) = if bits14 { (14, 10, 3, 512) } else { (12, 8, 4, 128) }; + let mut f = Fields::new(b); + let mut prev = [f.take(seed), f.take(seed)]; + put(out, 0, prev[0]); + put(out, 1, prev[1]); + for g in 0..groups { + let step = STEPS[f.take(2) as usize & 3]; + for k in 0..3 { + let i = 2 + 3 * g + k; + prev[i & 1] = next(prev[i & 1], f.take(code), step, bias); + put(out, i, prev[i & 1]); + } + } +} + +/// Decode a continuous stream of 16-byte blocks into `out`, in parallel by 0x4000-byte chunks (whole chunks +/// rotated when `rotated`). Blocks missing from `src` give zeros. +fn decode_blocks(src: &[u8], kind: Block, rotated: bool, out: &mut [u16]) { + let per = kind.per(); + out.par_chunks_mut(CHUNK / 16 * per).enumerate().for_each(|(c, o)| { + let raw = src.get(c * CHUNK..).unwrap_or_default(); + let raw = raw.get(..CHUNK).unwrap_or(raw); + let mut logical = Vec::new(); + let chunk = if rotated && raw.len() == CHUNK { + logical.extend_from_slice(raw.get(ROTATE..).unwrap_or_default()); + logical.extend_from_slice(raw.get(..ROTATE).unwrap_or_default()); + &logical + } else { + raw + }; + for (i, px) in o.chunks_mut(per).enumerate() { + match chunk.get(i * 16..i * 16 + 16).and_then(|s| <&[u8; 16]>::try_from(s).ok()) { + Some(b) => kind.decode(b, px), + None => px.fill(0), + } + } + }); +} + +fn decode_words(src: &[u8], shift: u32, out: &mut [u16]) { + out.par_chunks_mut(1 << 16).enumerate().for_each(|(k, o)| { + for (i, v) in o.iter_mut().enumerate() { + let at = ((k << 16) + i) * 2; + *v = src.get(at..at + 2).and_then(|b| <[u8; 2]>::try_from(b).ok()).map_or(0, |b| u16::from_le_bytes(b) >> shift); + } + }); +} + +/// A tag as a `u16` array: SHORT values, or the little-endian contents of an UNDEFINED / BYTE value. +fn u16_array(ifd: &Ifd, tag: u16, order: ByteOrder) -> Option> { + match ifd.value(tag)? { + Value::Undefined(b) | Value::Byte(b) => Some(b.as_chunks::<2>().0.iter().map(|c| order.u16(*c)).collect()), + Value::Short(v) => Some(v.clone()), + other => Some(other.to_u64_vec().into_iter().map(|v| v.min(u16::MAX as u64) as u16).collect()), + } +} + +/// The format-8 tile origins (tag `0x0045`): a count, then (x, y) pairs of `u16`. +fn origins(ifd: &Ifd, order: ByteOrder) -> Result> { + let a = u16_array(ifd, STRIP_LEFT, order).ok_or_else(|| RawError::Corrupt("RW2 format 8 without tag 0x0045".into()))?; + let n = a.first().map_or(0, |&n| n as usize); + let at = |i: usize| a.get(i).map(|&v| v as usize); + let values: Option> = (0..n).map(|i| Some((at(1 + 2 * i)?, at(2 + 2 * i)?))).collect(); + values.filter(|v| !v.is_empty() && v.len() <= MAX_STRIPS).ok_or_else(|| RawError::Corrupt("RW2 format 8: bad list in tag 0x0045".into())) +} + +/// A format-8 `u16` list (count, then the values; `wide`: 32-bit values as low/high pairs). +fn list(ifd: &Ifd, tag: u16, order: ByteOrder, wide: bool) -> Result> { + let a = u16_array(ifd, tag, order).ok_or_else(|| RawError::Corrupt(format!("RW2 format 8 without tag {tag:#06x}")))?; + let n = a.first().map_or(0, |&n| n as usize); + let at = |i: usize| a.get(i).map(|&v| v as u64); + let values: Option> = (0..n).map(|i| if wide { Some(at(1 + 2 * i)? | at(2 + 2 * i)? << 16) } else { at(1 + i) }).collect(); + values.filter(|v| !v.is_empty() && v.len() <= MAX_STRIPS).ok_or_else(|| RawError::Corrupt(format!("RW2 format 8: bad list in tag {tag:#06x}"))) +} + +/// One format-8 strip or tile: its bit stream and placement. +struct Strip<'a> { + data: &'a [u8], + left: usize, + top: usize, + width: usize, + height: usize, + bits: u64, +} + +/// The format-8 strips or tiles (together covering the whole image exactly once) and code table. +fn format8<'a>(ifd: &Ifd, order: ByteOrder, bytes: &'a [u8], w: usize, h: usize) -> Result<(Vec>, Huffman)> { + for (tag, want) in FORMAT8_FIXED { + if let Some(v) = u16_array(ifd, tag, order) + && v != want + { + return Err(RawError::Unsupported(format!("Panasonic raw format 8 with tag {tag:#06x} = {v:?}"))); + } + } + let (offsets, origins, bits) = (list(ifd, STRIP_OFFSETS, order, true)?, origins(ifd, order)?, list(ifd, STRIP_BITS, order, true)?); + let (widths, heights) = (list(ifd, STRIP_WIDTHS, order, false)?, list(ifd, STRIP_HEIGHTS, order, false)?); + let n = offsets.len(); + if [origins.len(), bits.len(), widths.len(), heights.len()].iter().any(|&l| l != n) { + return Err(RawError::Corrupt("RW2 format 8: strip lists of different lengths".into())); + } + let mut strips = Vec::with_capacity(n); + let mut area = 0usize; + for i in 0..n { + let ((left, top), width, height) = (origins[i], widths[i] as usize, heights[i] as usize); + let len = usize::try_from(bits[i].div_ceil(8)).unwrap_or(usize::MAX); + let data = usize::try_from(offsets[i]).ok().and_then(|o| bytes.get(o..o.checked_add(len)?)); + let inside = width > 0 && height > 0 && left + width <= w && top + height <= h; + let Some(data) = data.filter(|_| inside) else { + return Err(RawError::Corrupt(format!("RW2 format 8: strip {i} outside the image or file"))); + }; + // every sample costs at least one bit (no code is shorter), so a strip stating fewer bits than samples is + // corrupt: refusing it here keeps a few-byte file from sizing the full-sensor buffer + if bits[i] < (width * height) as u64 { + return Err(RawError::Corrupt(format!("RW2 format 8: strip {i} states {} bits for {} samples", bits[i], width * height))); + } + area += width * height; + strips.push(Strip { data, left, top, width, height, bits: bits[i] }); + } + // inside the image, no two overlapping and the areas adding up to the image: every sample has exactly one strip + let overlap = strips.iter().enumerate().any(|(i, a)| { + strips[i + 1..].iter().any(|b| a.left < b.left + b.width && b.left < a.left + a.width && a.top < b.top + b.height && b.top < a.top + a.height) + }); + if overlap || area != w * h { + return Err(RawError::Corrupt(format!("RW2 format 8: strips cover {area} of {} samples", w * h))); + } + // tag 0x0043 counts the rows of tiles; a file that says otherwise is not one we understand + let mut tops: Vec = strips.iter().map(|s| s.top).collect(); + tops.sort_unstable(); + tops.dedup(); + if let Some(v) = u16_array(ifd, TILE_ROWS, order) + && v != [tops.len() as u16] + { + return Err(RawError::Unsupported(format!("Panasonic raw format 8 with tag {TILE_ROWS:#06x} = {v:?} for {} tile rows", tops.len()))); + } + let codes = u16_array(ifd, HUFFMAN_CODES, order).unwrap_or_default(); + let pairs: Vec<(u32, u32)> = match (codes.first(), codes.get(1..35)) { + (Some(17), Some(c)) => c.as_chunks::<2>().0.iter().map(|&[l, c]| (l as u32, c as u32)).collect(), + _ => return Err(RawError::Corrupt("RW2 format 8: bad code table".into())), + }; + Ok((strips, code_table(&pairs)?)) +} + +/// The prefix-code table from (length, code) per difference category 0..=16. In 12- and 14-bit files the codes of the +/// categories the bit depth can't produce start with a shorter code (8 ones): read bit by bit, the shorter code wins +/// and the longer ones are unreachable. +fn code_table(pairs: &[(u32, u32)]) -> Result { + let mut table = Vec::with_capacity(pairs.len()); + for (i, &(la, ca)) in pairs.iter().enumerate() { + if la > 12 { + return Err(RawError::Unsupported(format!("RW2 format 8: {la}-bit code"))); + } + if la > 0 && ca >> la != 0 { + return Err(RawError::Corrupt("RW2 format 8: code longer than its length".into())); + } + let shadowed = pairs.iter().any(|&(lb, cb)| lb > 0 && lb < la && ca >> (la - lb) == cb); + let duplicate = pairs.iter().enumerate().any(|(j, &(lb, cb))| j != i && la > 0 && (lb, cb) == (la, ca)); + if duplicate { + return Err(RawError::Corrupt("RW2 format 8: two categories with the same code".into())); + } + table.push(if la == 0 || shadowed { (0, 0) } else { ((ca << (12 - la)) as u16, la as u8) }); + } + Huffman::new(&table) +} + +/// Decode one strip into its column range of each output row. +fn decode_strip(s: &Strip, codes: &Huffman, max: i32, mut rows: Vec<&mut [u16]>) -> Result<()> { + // the stream's bits are taken from each byte starting at bit 0: reversed, it is an ordinary MSB-first stream + let reversed: Vec = s.data.iter().map(|b| b.reverse_bits()).collect(); + let mut bits = Bits::new(&reversed); + let mut first = [0i32; 4]; + for pair in 0..s.height.div_ceil(2) { + let mut cell = first; + for cx in 0..s.width.div_ceil(2) { + for (k, v) in cell.iter_mut().enumerate() { + let d = diff(&mut bits, codes).ok_or_else(|| RawError::Corrupt(format!("RW2 format 8: invalid code in row {}", 2 * pair)))?; + *v = v.saturating_add(d); + // a value outside the sample depth is a corrupt stream, not a sample to clamp + let sample = u16::try_from(*v).ok().filter(|_| *v <= max); + let sample = + sample.ok_or_else(|| RawError::Corrupt(format!("RW2 format 8: sample {} outside the {max} range in row {}", *v, 2 * pair)))?; + if let Some(px) = rows.get_mut(2 * pair + k % 2).and_then(|r| r.get_mut(2 * cx + k / 2)) { + *px = sample; + } + } + if cx == 0 { + first = cell; + } + } + if bits.consumed_bits() as u64 > s.bits { + return Err(RawError::Corrupt(format!("RW2 format 8: strip data ends at row {} of {}", 2 * pair, s.height))); + } + } + Ok(()) +} + +/// `bits`: the sample depth; a decoded value above it (or below 0) makes the strip corrupt. +fn decode_strips(strips: &[Strip], codes: &Huffman, w: usize, bits: u32, out: &mut [u16]) -> Result<()> { + let max = (1i32 << bits.clamp(1, 16)) - 1; + let mut parts: Vec> = strips.iter().map(|s| Vec::with_capacity(s.height)).collect(); + for (y, row) in out.chunks_mut(w).enumerate() { + // the strips crossing this row, left to right, each take their columns off its front + let mut here: Vec = (0..strips.len()).filter(|&i| strips.get(i).is_some_and(|s| s.top <= y && y < s.top + s.height)).collect(); + here.sort_by_key(|&i| strips.get(i).map_or(0, |s| s.left)); + let (mut rest, mut x) = (row, 0usize); + for i in here { + let Some(s) = strips.get(i) else { continue }; + let skip = s.left.saturating_sub(x).min(rest.len()); + let (_, tail) = rest.split_at_mut(skip); + let (part, tail) = tail.split_at_mut(s.width.min(tail.len())); + if let Some(rows) = parts.get_mut(i) { + rows.push(part); + } + rest = tail; + x = s.left + s.width; + } + } + strips.par_iter().zip(parts).try_for_each(|(s, rows)| decode_strip(s, codes, max, rows)) +} + +/// Replace samples of 0 inside `area` (mapped-out defects) by the mean of the non-zero samples of the same colour two +/// pixels away. +fn fill_defects(data: &mut [u16], w: usize, area: Rect) { + for y in area.y..area.y + area.height { + let Some(row) = data.get(y * w + area.x..y * w + area.x + area.width) else { break }; + if !row.contains(&0) { + continue; + } + for x in area.x..area.x + area.width { + let i = y * w + x; + if data.get(i) != Some(&0) { + continue; + } + let near = [(x.checked_sub(2), Some(y)), (Some(x + 2).filter(|&x| x < w), Some(y)), (Some(x), y.checked_sub(2)), (Some(x), Some(y + 2))]; + let (mut sum, mut n) = (0u32, 0u32); + for (nx, ny) in near { + if let (Some(nx), Some(ny)) = (nx, ny) + && let Some(&v) = data.get(ny * w + nx).filter(|&&v| v != 0) + { + sum += v as u32; + n += 1; + } + } + if n > 0 + && let Some(v) = data.get_mut(i) + { + *v = (sum / n) as u16; + } + } + } +} + +fn cfa(code: u16) -> Cfa { + Cfa::bayer_static(match code { + 2 => "GRBG", + 3 => "GBRG", + 4 => "BGGR", + _ => "RGGB", + }) +} + +pub(crate) fn decode(bytes: &[u8], mode: Mode, limits: &Limits) -> Result { + let tiff = Tiff::parse(bytes)?; + let ifd0 = tiff.ifds.first().ok_or_else(|| RawError::Corrupt("RW2 without IFD0".into()))?; + let get = |tag: u16| ifd0.u64(tag).and_then(|v| usize::try_from(v).ok()); + let (w, h) = (get(SENSOR_WIDTH).unwrap_or(0), get(SENSOR_HEIGHT).unwrap_or(0)); + if w == 0 || h == 0 { + return Err(RawError::Corrupt("RW2 without sensor size".into())); + } + let n = w.checked_mul(h).filter(|n| *n <= limits.samples()).ok_or(RawError::Limit("image too large"))?; + let bits = ifd0.u32(BITS).unwrap_or(12); + let encoding = encoding(ifd0, bits)?; + let src = || { + let off = get(RAW_OFFSET).or_else(|| get(t::STRIP_OFFSETS)).ok_or_else(|| RawError::Corrupt("RW2 without raw data offset".into()))?; + bytes.get(off..).ok_or_else(|| RawError::Corrupt("RW2 raw data outside file".into())) + }; + let truncated = |have: usize, need: usize| RawError::Corrupt(format!("RW2 raw data truncated ({have} of {need} bytes)")); + let mut data = Vec::new(); + match encoding { + Encoding::Words => { + let src = src()?; + let need = n * 2; + if src.len() < need { + return Err(truncated(src.len(), need)); + } + if mode == Mode::Full { + data = vec![0u16; n]; + decode_words(src, 16u32.saturating_sub(bits), &mut data); + } + } + Encoding::Blocks { kind, rotated } => { + let src = src()?; + let need = n.div_ceil(kind.per()) * 16; + if src.len() < need { + return Err(truncated(src.len(), need)); + } + if mode == Mode::Full { + data = vec![0u16; n]; + decode_blocks(src, kind, rotated, &mut data); + } + } + Encoding::Strips => { + let (strips, codes) = format8(ifd0, tiff.order, bytes, w, h)?; + if mode == Mode::Full { + data = vec![0u16; n]; + decode_strips(&strips, &codes, w, bits, &mut data)?; + } + } + } + + let active = match (get(BORDER_TOP), get(BORDER_LEFT), get(BORDER_BOTTOM), get(BORDER_RIGHT)) { + (Some(tp), Some(l), Some(b), Some(r)) if b > tp && r > l => Rect::new(l, tp, r - l, b - tp).clipped(w, h), + _ => Rect::new(0, 0, w, h), + }; + let active = if active.width == 0 || active.height == 0 { Rect::new(0, 0, w, h) } else { active }; + let crop = match CROP.map(get) { + [Some(tp), Some(l), Some(b), Some(r)] + if b > tp && r > l && l >= active.x && tp >= active.y && r <= active.x + active.width && b <= active.y + active.height => + { + Rect::new(l - active.x, tp - active.y, r - l, b - tp) + } + _ => Rect::new(0, 0, active.width, active.height), + }; + if !data.is_empty() { + fill_defects(&mut data, w, active); + } + + let cfa = cfa(ifd0.u16(CFA_PATTERN).unwrap_or(1)); + let tags = BLACK_RGB.map(|tag| ifd0.f64(tag).filter(|v| v.is_finite() && *v >= 0.0).unwrap_or(0.0) as f32); + let offset = match encoding { + Encoding::Blocks { kind: Block::Format4 | Block::Format6 { .. }, .. } => CODEC_BLACK_OFFSET, + _ => 0.0, + }; + let rgb = if tags.iter().all(|&v| v == 0.0) { [DEFAULT_BLACK; 3] } else { tags.map(|v| v + offset) }; + // per CFA position, anchored at the active area + let values = cfa.shifted(active.x, active.y).pattern.iter().map(|&c| rgb[c as usize]).collect(); + let black = BlackLevel { repeat_rows: 2, repeat_cols: 2, values, ..Default::default() }; + let full = ((1u32 << bits.clamp(1, 16)) - 1) as f32; + let white = if data.is_empty() { + // headers only (white isn't part of `RawInfo`): the recorded limit as a placeholder + LINEARITY_LIMIT.iter().filter_map(|&tag| ifd0.f64(tag)).map(|v| v as f32).filter(|v| *v > 0.0 && *v <= full).fold(full, f32::min) + } else { + white_from_data(&active_rows(&data, w, active), bits) + }; + let white = if white > black.mean() { white } else { full }; + let wb = match WB_RGB.map(|tag| ifd0.f64(tag)) { + [Some(r), Some(g), Some(b)] if r > 0.0 && g > 0.0 && b > 0.0 => Some([(r / g) as f32, 1.0, (b / g) as f32]), + _ => None, + }; + let img = RawImage { + format: RawFormat::Rw2, + width: w, + height: h, + cpp: 1, + data: RawData::U16(data), + cfa: Some(cfa), + bits, + black, + white: vec![white], + active_area: active, + crop, + orientation: Orientation::from_exif(ifd0.u16(t::ORIENTATION).unwrap_or(1)), + color: ColorData::default(), + wb_multipliers: wb, + linearized: false, + pixel_aspect: 1.0, + opcodes: OpcodeLists { list3: distortion(ifd0, active).into_iter().collect(), ..Default::default() }, + camera: tiff.ifds.first().map(Camera::from_ifd0).unwrap_or_default(), + cinema: Default::default(), + }; + img.validate_for(mode)?; + Ok(img) +} + +/// The camera's distortion correction (tag `0x0119`, see the module docs) as the equivalent DNG `WarpRectilinear` +/// opcode for `OpcodeList3` (relative to the active area), so it reaches the same "profile corrections" as a DNG's +/// embedded warp. `None` when the camera recorded the correction as off, when it is the identity, or when the +/// values are implausible (the photo is then shown uncorrected rather than wrongly warped). +fn distortion(ifd: &Ifd, active: Rect) -> Option { + let b = ifd.bytes(DISTORTION_INFO)?; + let word = |i: usize| -> Option<[u8; 2]> { b.get(2 * i..2 * i + 2)?.try_into().ok() }; + let signed = |i: usize| word(i).map(|w| i16::from_le_bytes(w) as f64); + if u16::from_le_bytes(word(7)?) & 0x000f != 1 { + return None; + } + let (k2, k4, k6) = (signed(8)? / DISTORTION_ONE, signed(4)? / DISTORTION_ONE, signed(11)? / DISTORTION_ONE); + let zoom = 1.0 + signed(5)? / DISTORTION_ONE; + let radius = u16::from_le_bytes(word(12)?) as f64; + if k2 == 0.0 && k4 == 0.0 && k6 == 0.0 && zoom == 1.0 { + return None; + } + // DNG's normalisation: the distance from the centre to the farthest corner, in pixel-index units + let half_diag = (active.width as f64 - 1.0).hypot(active.height as f64 - 1.0) / 2.0; + if !(half_diag >= 1.0 && (0.5 * half_diag..=2.0 * half_diag).contains(&radius) && (0.75..=1.25).contains(&zoom)) { + return None; + } + // the camera's correction: a source point at radius ρ (in units of `radius`) moves to ρ·g(ρ) + let corrected = |r: f64| { + let r2 = r * r; + r * (1.0 + r2 * (k2 + r2 * (k4 + r2 * k6))) + }; + // work in DNG-normalised units (corner distance = 1), where the output corners sit at radius `zoom` + let scale = radius / half_diag; + let k = crate::shared::warp_from_forward_map(&|s| corrected(s / scale) * scale, zoom, DISTORTION_MAX_FIT_ERROR)?; + Some(Opcode::WarpRectilinear { planes: vec![[k[0], k[1], k[2], k[3], 0.0, 0.0]], center: [0.5, 0.5] }) +} + +#[cfg(test)] +mod tests { + use super::*; + use craft_tiff::{IfdBuilder, TiffWriter}; + + /// Writes bit fields into a 16-byte block from bit 127 down (the inverse of [`Fields`]). + struct FieldWriter { + word: u128, + left: u32, + } + + impl FieldWriter { + fn new() -> Self { + FieldWriter { word: 0, left: 128 } + } + fn put(&mut self, n: u32, v: u32) -> &mut Self { + assert!(v < 1 << n, "{v} doesn't fit {n} bits"); + self.left -= n; + self.word |= (v as u128) << self.left; + self + } + fn block(&self) -> [u8; 16] { + self.word.to_le_bytes() + } + } + + /// A seed of format 4: 8 bits `h`, 4 bits `l`. + fn seed(f: &mut FieldWriter, v: u32) { + f.put(8, v >> 4).put(4, v & 15); + } + + fn decode4(b: [u8; 16]) -> [u16; 14] { + let mut out = [0u16; 14]; + format4(&b, &mut out); + out + } + + #[test] + fn format4_known_blocks() { + // every branch: differences at each scale, a zero code, absolute scale-3 codes keeping the low bits + let mut f = FieldWriter::new(); + seed(&mut f, 1000); + seed(&mut f, 2000); + f.put(2, 0).put(8, 130).put(8, 120).put(8, 0); // 1002, 1992, 1002 (kept) + f.put(2, 1).put(8, 129).put(8, 100).put(8, 128); // 1994, 946, 1994 + f.put(2, 2).put(8, 200).put(8, 64).put(8, 131); // 1234, 1738, 1246 + f.put(2, 3).put(8, 50).put(8, 255).put(8, 0); // 50·16 + 1738 % 16 = 810, 4080 + 14 = 4094, 810 (kept) + assert_eq!(f.left, 0); + assert_eq!(decode4(f.block()), [1000, 2000, 1002, 1992, 1002, 1994, 946, 1994, 1234, 1738, 1246, 810, 4094, 810]); + + // a deferred seed (defect marker, 8 bits only) and absolute codes because the previous value is too small + let mut f = FieldWriter::new(); + f.put(8, 0); // pixel 0: defect; the next even pixel brings the seed + seed(&mut f, 325); + f.put(2, 1); + seed(&mut f, 481); // pixel 2: the even seed + f.put(8, 140).put(8, 10); // 325 + 24 = 349, 481 − 236 = 245 + f.put(2, 2).put(8, 100).put(8, 129).put(8, 129); // < 512: 400 + 349 % 4 = 401, 516 + 245 % 4 = 517, 517 + f.put(2, 0).put(8, 127).put(8, 0).put(8, 255); // 516, 517 (kept), 643 + f.put(2, 0).put(8, 1).put(8, 128).put(8, 128); // 390, 643, 390 + assert_eq!(f.left, 0); + assert_eq!(decode4(f.block()), [0, 325, 481, 349, 245, 401, 517, 517, 516, 517, 643, 390, 643, 390]); + + // a colour without a seed until pixel 12, which then always reads its 4 low bits + let mut f = FieldWriter::new(); + f.put(8, 0); + seed(&mut f, 1600); + for g in 0..4 { + f.put(2, 0); + for k in 0..3 { + let i = 2 + 3 * g + k; + if i % 2 == 1 { + f.put(8, 128); + } else if i < 12 { + f.put(8, 0); + } else { + f.put(8, 0).put(4, 7); + } + } + } + assert_eq!(f.left, 0); + assert_eq!(decode4(f.block()), [0, 1600, 0, 1600, 0, 1600, 0, 1600, 0, 1600, 0, 1600, 7, 1600]); + } + + #[test] + fn format6_known_blocks() { + let mut f = FieldWriter::new(); + f.put(14, 5000).put(14, 6000); + f.put(2, 0).put(10, 520).put(10, 500).put(10, 512); // 5008, 5988, 5008 + f.put(2, 3).put(10, 1023).put(10, 300).put(10, 0); // 16368 + 5988 % 16, 4800 + 5008 % 16, kept + f.put(2, 1).put(10, 600).put(10, 400).put(10, 10); // 4976, 16148, 3972 + let mut out = [0u16; 11]; + format6(&f.block(), true, &mut out); + assert_eq!(out, [5000, 6000, 5008, 5988, 5008, 16372, 4800, 16372, 4976, 16148, 3972]); + + let mut f = FieldWriter::new(); + f.put(12, 300).put(12, 4000); + f.put(2, 2).put(8, 140).put(8, 100).put(8, 130); // 300 < 512: 560 + 0; 3888; 568 + f.put(2, 0).put(8, 129).put(8, 127).put(8, 128); // 3889, 567, 3889 + f.put(2, 3).put(8, 10).put(8, 255).put(8, 20); // 160 + 7, 4080 + 1, 320 + 7 + f.put(2, 1).put(8, 128).put(8, 130).put(8, 0); // 4081, 331, 4081 (kept) + assert_eq!(f.left, 0); + let mut out = [0u16; 14]; + format6(&f.block(), false, &mut out); + assert_eq!(out, [300, 4000, 560, 3888, 568, 3889, 567, 3889, 167, 4081, 327, 4081, 331, 4081]); + } + + /// A camera-like format 4/6 encoder: per group the smallest scale that codes the targets exactly, else scale 3. + /// Returns the block and what a decoder must reconstruct. + fn encode_block(target: &[u16], fmt6_bits: Option) -> ([u8; 16], Vec) { + let (seed_bits, code_bits, bias) = match fmt6_bits { + Some(14) => (14, 10, 512), + Some(_) => (12, 8, 128), + None => (12, 8, 128), + }; + let n = target.len(); + let mut f = FieldWriter::new(); + let mut prev = [target[0] as i32, target[1] as i32]; + let mut out = vec![target[0], target[1]]; + for &s in &target[..2] { + if fmt6_bits.is_some() { + f.put(seed_bits, s as u32); + } else { + seed(&mut f, s as u32); + } + } + for g in 0..(n - 2) / 3 { + let idx: Vec = (0..3).map(|k| 2 + 3 * g + k).collect(); + let fits = |s: usize| -> Option> { + let (m, absolute) = STEPS[s]; + let mut p = prev; + let mut codes = vec![]; + for &i in &idx { + let (v, c) = (target[i] as i32, i & 1); + let j = if !absolute && p[c] >= bias * m { + let d = v - p[c]; + (d % m == 0 && (d / m).abs() < bias).then_some(d / m + bias)? + } else { + // nearest step, as a camera would + let j = (v - (p[c] & (m - 1)) + m / 2) / m; + (absolute || (v - (p[c] & (m - 1))) % m == 0).then_some(j.clamp(1, (1 << code_bits) - 1))? + }; + p[c] = next(p[c], j, (m, absolute), bias); + codes.push(j); + } + Some(codes) + }; + let s = (0..4).find(|&s| fits(s).is_some()).unwrap_or(3); + let codes = fits(s).unwrap(); + f.put(2, s as u32); + for (k, &i) in idx.iter().enumerate() { + f.put(code_bits, codes[k] as u32); + prev[i & 1] = next(prev[i & 1], codes[k], STEPS[s], bias); + out.push(prev[i & 1] as u16); + } + } + (f.block(), out) + } + + /// A smooth synthetic Bayer-like image with some sharp edges, values `lo..hi`. + fn image(w: usize, h: usize, lo: u32, hi: u32) -> Vec { + (0..w * h) + .map(|i| { + let (x, y) = (i % w, i / w); + let smooth = lo + (hi - lo) * ((x * 7 + y * 3) % 997) as u32 / 997 / 2; + let edge = if (x / 37 + y / 23) % 3 == 0 { (hi - lo) / 3 } else { 0 }; + (smooth + edge + (x % 2 * 40) as u32).min(hi) as u16 + }) + .collect() + } + + /// Rotate every whole 0x4000-byte chunk the way formats 4 and 5 store them. + fn rotate(mut stream: Vec) -> Vec { + stream.resize(stream.len().div_ceil(CHUNK) * CHUNK, 0); + stream.chunks(CHUNK).flat_map(|c| c[CHUNK - ROTATE..].iter().chain(&c[..CHUNK - ROTATE]).copied().collect::>()).collect() + } + + /// Pack samples LSB-first into 16-byte blocks (formats 2, 5 and 7). + fn pack(px: &[u16], bits: u32) -> Vec { + let per = (128 / bits) as usize; + px.chunks(per).flat_map(|blk| blk.iter().enumerate().fold(0u128, |w, (i, &v)| w | (v as u128) << (i as u32 * bits)).to_le_bytes()).collect() + } + + /// A minimal RW2: IFD0 with `tags` plus sensor size, borders (2 rows, 4 columns at the top/left) and Make, the + /// payload at the end. `at_payload` adds tags that hold the payload's offset. + fn file(w: u16, h: u16, tags: &[(u16, Value)], payload: &[u8], at_payload: &dyn Fn(u32) -> Vec<(u16, Value)>) -> Vec { + let build = |off: u32| { + let mut ifd = IfdBuilder::new(); + for (tag, v) in [(SENSOR_WIDTH, w), (SENSOR_HEIGHT, h), (BORDER_TOP, 2), (BORDER_LEFT, 4), (BORDER_BOTTOM, h), (BORDER_RIGHT, w)] { + ifd.set(tag, Value::Short(vec![v])); + } + ifd.set(t::MAKE, Value::Ascii("Panasonic".into())); + for (tag, v) in tags.iter().cloned().chain(at_payload(off)) { + ifd.set(tag, v); + } + let mut b = TiffWriter::new(ByteOrder::Little, false).write(&[ifd]).unwrap(); + b[2] = 0x55; + b + }; + let off = build(0).len().div_ceil(16) * 16; + let mut b = build(off as u32); + b.resize(off, 0); + b.extend_from_slice(payload); + b + } + + fn short(tag: u16, v: u16) -> (u16, Value) { + (tag, Value::Short(vec![v])) + } + + /// Tags of a format-`format`, `bits`-bit file with black levels `black` and the payload at `RawDataOffset`. + fn rw2(w: u16, h: u16, bits: u16, format: u16, black: u16, payload: &[u8]) -> Vec { + let mut tags = vec![short(CFA_PATTERN, 1), short(BITS, bits), short(RAW_FORMAT, format)]; + tags.extend(BLACK_RGB.iter().enumerate().map(|(i, &tag)| short(tag, black + i as u16))); + tags.extend(WB_RGB.iter().zip([512, 256, 384]).map(|(&tag, v)| short(tag, v))); + file(w, h, &tags, payload, &|off| vec![(RAW_OFFSET, Value::Long(vec![off]))]) + } + + fn samples(r: &RawImage) -> &[u16] { + match &r.data { + RawData::U16(d) => d, + RawData::F32(_) => panic!("float samples"), + } + } + + #[test] + fn format4_and_6_round_trip_through_files() { + // format 4: 14 pixels per block, rotated chunks, blocks continuous across rows (here 70 = 5 blocks per row) + let (w, h) = (70usize, 300usize); + let target = image(w, h, 200, 4000); + let (mut stream, mut want) = (vec![], vec![]); + for px in target.chunks(14) { + let (b, out) = encode_block(px, None); + stream.extend_from_slice(&b); + want.extend(out); + } + assert!(want.iter().zip(&target).all(|(a, b)| a.abs_diff(*b) <= 8), "the encoder is too lossy for the test image"); + let bytes = rw2(w as u16, h as u16, 12, 4, 128, &rotate(stream)); + let r = crate::decode(&bytes).unwrap(); + assert_eq!(samples(&r), want); + assert_eq!(crate::probe_info(&bytes).unwrap(), r.info()); + assert_eq!(r.black.values, vec![143.0; 4].into_iter().zip([0., 1., 1., 2.]).map(|(b, o)| b + o).collect::>()); + + // format 6: not rotated; 11 (14-bit) or 14 (12-bit) pixels per block + for (bits, per, top) in [(14u16, 11usize, 16000u32), (12, 14, 4000)] { + let w = per * 6; + let target = image(w, 40, 600, top); + let (mut stream, mut want) = (vec![], vec![]); + for px in target.chunks(per) { + let (b, out) = encode_block(px, Some(bits as u32)); + stream.extend_from_slice(&b); + want.extend(out); + } + let bytes = rw2(w as u16, 40, bits, 6, 512, &stream); + let r = crate::decode(&bytes).unwrap_or_else(|e| panic!("{bits}: {e}")); + assert_eq!(samples(&r), want, "{bits}-bit"); + assert_eq!(r.black.at(0, 0, 0, 1), 527.0, "format 6 stores samples 15 above the black level"); + } + } + + /// A tag-0x0119 block: 16 little-endian words, the given ones set. + fn distortion_block(words: &[(usize, i16)]) -> Value { + let mut w = [0i16; 16]; + for &(i, v) in words { + w[i] = v; + } + Value::Undefined(w.iter().flat_map(|v| v.to_le_bytes()).collect()) + } + + fn distorted(info: Option) -> Vec { + let (w, h) = (90usize, 60usize); + let px: Vec = (0..w * h).map(|i| (i % 4000) as u16 + 1).collect(); + let mut tags = vec![short(CFA_PATTERN, 1), short(BITS, 12), short(RAW_FORMAT, 5)]; + tags.extend(info.map(|v| (DISTORTION_INFO, v))); + let bytes = file(w as u16, h as u16, &tags, &rotate(pack(&px, 12)), &|off| vec![(RAW_OFFSET, Value::Long(vec![off]))]); + let r = crate::decode(&bytes).unwrap(); + assert_eq!(crate::probe_info(&bytes).unwrap().opcodes, r.opcodes, "headers-only probe reads the same correction"); + r.opcodes.list3 + } + + /// Known answer: correction k2 = 4096/32768 = 1/8 (word 8), zoom 1 − 1024/32768 (word 5), radius 51 px (word 12) + /// on an 86 × 58 active area. A sensor point at ρ appears at ρ + ρ³/8, so the output corner (corrected radius + /// t = half diagonal · zoom / 51) shows the root of ρ³ + 8ρ − 8t = 0, here by Cardano's formula. + #[test] + fn distortion_correction_becomes_a_rectilinear_warp() { + let on = |more: &[(usize, i16)]| { + let mut w = vec![(7, 0x0001), (8, 4096), (5, -1024), (12, 51)]; + w.extend_from_slice(more); + distortion_block(&w) + }; + let ops = distorted(Some(on(&[]))); + let [Opcode::WarpRectilinear { planes, center }] = ops.as_slice() else { panic!("{ops:?}") }; + assert_eq!((planes.len(), *center), (1, [0.5, 0.5])); + let k = planes[0]; + assert_eq!((k[4], k[5]), (0.0, 0.0), "radial only"); + let half_diag = 85f64.hypot(57.0) / 2.0; + for u in [0.1, 0.25, 0.5, 0.75, 1.0] { + let t = u * half_diag * (1.0 - 1024.0 / 32768.0) / 51.0; + let (p, q) = (8.0, -8.0 * t); + let d = (q * q / 4.0 + p * p * p / 27.0).sqrt(); + let rho = (-q / 2.0 + d).cbrt() + (-q / 2.0 - d).cbrt(); + let want = rho * 51.0 / half_diag; + let got = u * (k[0] + k[1] * u * u + k[2] * u.powi(4) + k[3] * u.powi(6)); + // the 4-term odd polynomial approximates the exact inverse (5e-5 of the half diagonal ≈ 0.2 px at 20 MP) + assert!((got - want).abs() < 5e-5, "u {u}: source radius {got}, want {want}"); + } + + // no correction: no block, recorded as off, or the identity + assert!(distorted(None).is_empty()); + assert!(distorted(Some(distortion_block(&[(7, 0x0000), (8, 4096), (12, 51)]))).is_empty(), "off"); + assert!(distorted(Some(distortion_block(&[(7, 0x7f00u16 as i16), (8, 4096), (12, 51)]))).is_empty(), "off, other bits set"); + assert!(distorted(Some(distortion_block(&[(7, 0x0001), (12, 51)]))).is_empty(), "identity"); + // hostile values are ignored, never a panic or a wild warp + for (what, block) in [ + ("truncated", Value::Undefined(vec![1, 0, 2, 0, 3])), + ("radius 0", on(&[(12, 0)])), + ("radius far beyond the sensor", on(&[(12, -1)])), + ("zoom 0", on(&[(5, -32768)])), + ("folds over", on(&[(8, -32768)])), + ("folds over at high order", on(&[(11, -32768), (4, -32768)])), + ("runs away", on(&[(8, 32767), (4, 32767), (11, 32767)])), + ] { + let ops = distorted(Some(block)); + assert!(ops.iter().all(|op| !matches!(op, Opcode::WarpRectilinear { .. })) || what == "runs away", "{what}: {ops:?}"); + if let [Opcode::WarpRectilinear { planes, .. }] = ops.as_slice() { + assert!(planes[0].iter().all(|c| c.is_finite()), "{what}"); + } + } + } + + #[test] + fn packed_and_word_formats() { + let (w, h) = (90usize, 200usize); // > one chunk of data + for (format, bits, rotated) in [(5u16, 12u32, true), (5, 14, true), (7, 14, false), (2, 12, false)] { + let px: Vec = (0..w * h).map(|i| (((i * 2654435761usize) >> 7) as u16 & ((1 << bits) - 1)).max(1)).collect(); + let stream = pack(&px, bits); + let bytes = rw2(w as u16, h as u16, bits as u16, format, 100, &if rotated { rotate(stream) } else { stream }); + assert_eq!(crate::probe(&bytes), Some(RawFormat::Rw2)); + let r = crate::decode(&bytes).unwrap_or_else(|e| panic!("format {format}, {bits}-bit: {e}")); + assert_eq!(crate::probe_info(&bytes).unwrap(), r.info()); + assert_eq!(samples(&r), px, "format {format}, {bits}-bit"); + assert_eq!(r.active_area, Rect::new(4, 2, w - 4, h - 2)); + assert_eq!(r.cfa.as_ref().unwrap().name(), "RGGB"); + // packed formats store samples at the recorded black level, per colour, anchored at the active area + assert_eq!([r.black.at(0, 0, 0, 1), r.black.at(1, 0, 0, 1), r.black.at(1, 1, 0, 1)], [100.0, 101.0, 102.0]); + assert_eq!(r.wb_multipliers, Some([2.0, 1.0, 1.5])); + } + // the oldest files: 16-bit words with the sample left-justified, no black level recorded + for (compression, bits) in [(34828u16, 12u32), (34830, 14)] { + let px: Vec = (0..w * h).map(|i| 100 + (i * 7 % 3000) as u16).collect(); + let payload: Vec = px.iter().flat_map(|&v| (v << (16 - bits)).to_le_bytes()).collect(); + let tags = [short(BITS, bits as u16), short(COMPRESSION, compression), short(CFA_PATTERN, 4)]; + let bytes = file(w as u16, h as u16, &tags, &payload, &|off| vec![(t::STRIP_OFFSETS, Value::Long(vec![off]))]); + let r = crate::decode(&bytes).unwrap_or_else(|e| panic!("{compression}: {e}")); + assert_eq!(samples(&r), px); + assert_eq!(r.black.mean(), DEFAULT_BLACK); + assert_eq!(r.white_at(0), ((1 << bits) - 1) as f32, "no saturation plateau: the full range"); + } + } + + #[test] + fn levels_crop_defects_and_cfa() { + let (w, h) = (90usize, 60usize); + let mut px: Vec = (0..w * h).map(|i| 300 + (i % 7) as u16).collect(); + // a saturation plateau sets the white level; zeros inside the active area are defects + px[10 * w..20 * w].fill(3900); + px[30 * w + 40] = 0; + px[1] = 0; // in the border: left alone + let mut tags = vec![short(CFA_PATTERN, 2), short(BITS, 12), short(RAW_FORMAT, 5)]; + tags.extend([0x002f, 0x0030, 0x0031, 0x0032].iter().zip([12u16, 10, 52, 80]).map(|(&tag, v)| short(tag, v))); + tags.extend(LINEARITY_LIMIT.map(|tag| short(tag, 2000))); // recorded limits are not used + let bytes = file(w as u16, h as u16, &tags, &rotate(pack(&px, 12)), &|off| vec![(RAW_OFFSET, Value::Long(vec![off]))]); + let r = crate::decode(&bytes).unwrap(); + let d = samples(&r); + assert_eq!(d[30 * w + 40], ((px[30 * w + 38] as u32 + px[30 * w + 42] as u32 + px[28 * w + 40] as u32 + px[32 * w + 40] as u32) / 4) as u16); + assert_eq!(d[1], 0); + assert!((3800.0..=3900.0).contains(&r.white_at(0)), "{}", r.white_at(0)); + assert_eq!(r.black.mean(), DEFAULT_BLACK, "no black recorded"); + assert_eq!(r.crop, Rect::new(6, 10, 70, 40), "crop tags are sensor coordinates; the crop is relative to the active area"); + assert_eq!(r.info().developed_size(), (70, 40)); + // the CFA tag describes sample (0, 0) of the data, not the active area's origin + assert_eq!(r.cfa.as_ref().unwrap().name(), "GRBG"); + let mut tags = vec![short(CFA_PATTERN, 2), short(BITS, 12), short(RAW_FORMAT, 5)]; + tags.extend([0x002f, 0x0030, 0x0031, 0x0032].map(|tag| short(tag, 0))); // some bodies write zeros + let bytes = file(w as u16, h as u16, &tags, &rotate(pack(&px, 12)), &|off| vec![(RAW_OFFSET, Value::Long(vec![off]))]); + assert_eq!(crate::decode(&bytes).unwrap().crop, Rect::new(0, 0, w - 4, h - 2)); + } + + /// Format 8 code tables as two cameras write them (16-bit DC-GH6; 14-bit DC-S5M2, whose codes for categories 15 + /// and 16 start with the 8-bit code of category 14 and are unreachable). + const TABLE16: [(u32, u32); 17] = [ + (10, 1022), + (11, 2046), + (8, 254), + (9, 510), + (7, 126), + (4, 14), + (4, 12), + (3, 4), + (3, 2), + (2, 0), + (3, 3), + (3, 5), + (4, 13), + (5, 30), + (6, 62), + (12, 4094), + (12, 4095), + ]; + type Tags = Vec<(u16, Value)>; + + const TABLE14: [(u32, u32); 17] = [ + (6, 62), + (7, 126), + (6, 61), + (5, 28), + (4, 12), + (3, 4), + (3, 2), + (2, 0), + (3, 3), + (3, 5), + (4, 13), + (5, 29), + (6, 60), + (8, 254), + (8, 255), + (12, 4094), + (12, 4095), + ]; + + /// Encode a strip (the inverse of [`decode_strip`]): returns the stored bytes and the number of bits. + fn encode_strip(px: &[u16], w: usize, h: usize, table: &[(u32, u32); 17], padding_rows: usize) -> (Vec, u64) { + let mut bits: Vec = vec![]; + let mut put = |n: u32, v: u32| bits.extend((0..n).rev().map(|i| v >> i & 1 == 1)); + let mut first = [0i32; 4]; + let at = |x: usize, y: usize| if y < h { px[y * w + x] as i32 } else { 1000 }; + for pair in 0..(h / 2 + padding_rows) { + let mut cell = first; + for cx in 0..w / 2 { + for (k, p) in cell.iter_mut().enumerate() { + let v = at(2 * cx + k / 2, 2 * pair + k % 2); + let d = v - *p; + let c = 32 - d.unsigned_abs().leading_zeros(); + let (len, code) = table[c as usize]; + put(len, code); + put(c, if d >= 0 { d as u32 } else { (d + (1 << c) - 1) as u32 }); + *p = v; + } + if cx == 0 { + first = cell; + } + } + } + let n = bits.len() as u64; + // stored with each byte's bits taken from bit 0 up + (bits.chunks(8).map(|b| b.iter().enumerate().fold(0u8, |a, (i, &x)| a | (x as u8) << i)).collect(), n) + } + + /// A format-8 file with the given strip widths; `tweak` edits the tags before writing. + fn format8_file(w: usize, h: usize, px: &[u16], widths: &[usize], table: &[(u32, u32); 17], tweak: &dyn Fn(&mut Tags)) -> Vec { + format8_tiled_file(w, px, widths, &[h], table, tweak) + } + + /// A format-8 file tiled in rows of the given heights, each cut into the given widths. + fn format8_tiled_file(w: usize, px: &[u16], widths: &[usize], rows: &[usize], table: &[(u32, u32); 17], tweak: &dyn Fn(&mut Tags)) -> Vec { + let h: usize = rows.iter().sum(); + let mut payload = vec![]; + let mut strips = vec![]; // (offset in payload, bits, left, width, top, height) + let mut top = 0; + for &th in rows { + let mut left = 0; + for &sw in widths { + let sub: Vec = (top..top + th).flat_map(|y| px[y * w + left..y * w + left + sw].to_vec()).collect(); + let (data, bits) = encode_strip(&sub, sw, th, table, 2); + strips.push((payload.len(), bits, left, sw, top, th)); + payload.extend(data); + payload.resize(payload.len().div_ceil(16) * 16, 0); + left += sw; + } + top += th; + } + let u16s = |v: Vec| undefined(&v); + let wide = |v: Vec| u16s(std::iter::once(v.len() as u16).chain(v.iter().flat_map(|&x| [x as u16, (x >> 16) as u16])).collect()); + let narrow = |v: Vec| u16s(std::iter::once(v.len() as u16).chain(v.iter().map(|&x| x as u16)).collect()); + let mut tags = vec![short(CFA_PATTERN, 1), short(BITS, 14), short(RAW_FORMAT, 8)]; + tags.extend(BLACK_RGB.map(|tag| short(tag, 512))); + tags.push((HUFFMAN_CODES, u16s(std::iter::once(17).chain(table.iter().flat_map(|&(l, c)| [l as u16, c as u16])).collect()))); + for (tag, v) in FORMAT8_FIXED { + tags.push(if v.len() == 1 { short(tag, v[0]) } else { (tag, u16s(v.to_vec())) }); + } + tags.push(short(TILE_ROWS, rows.len() as u16)); + tags.push((STRIP_LEFT, u16s(std::iter::once(strips.len() as u16).chain(strips.iter().flat_map(|s| [s.2 as u16, s.4 as u16])).collect()))); + tags.push((STRIP_BITS, wide(strips.iter().map(|s| s.1).collect()))); + tags.push((STRIP_WIDTHS, narrow(strips.iter().map(|s| s.3 as u64).collect()))); + tags.push((STRIP_HEIGHTS, narrow(strips.iter().map(|s| s.5 as u64).collect()))); + tweak(&mut tags); + let offsets: Vec = strips.iter().map(|s| s.0).collect(); + file(w as u16, h as u16, &tags, &payload, &|off| { + vec![(RAW_OFFSET, Value::Long(vec![off])), (STRIP_OFFSETS, wide(offsets.iter().map(|&o| off as u64 + o as u64).collect()))] + }) + } + + #[test] + fn format8_strips_round_trip() { + let (w, h) = (64usize, 30usize); + let px = image(w, h, 520, 16000); + for (table, widths) in [(&TABLE14, vec![64]), (&TABLE14, vec![30, 34]), (&TABLE16, vec![20, 22, 22])] { + let bytes = format8_file(w, h, &px, &widths, table, &|_| {}); + let r = crate::decode(&bytes).unwrap_or_else(|e| panic!("{widths:?}: {e}")); + assert_eq!(samples(&r), px, "{widths:?}"); + assert_eq!(crate::probe_info(&bytes).unwrap(), r.info()); + assert_eq!(r.black.mean(), 512.0, "format 8 stores samples at the recorded black level"); + } + } + + /// A format-8 sample above the file's bit depth is a corrupt stream: refused, not clamped or kept. The top of + /// the range itself decodes. + #[test] + fn format8_samples_above_the_bit_depth_are_refused() { + let (w, h) = (64usize, 30usize); + for depth in [12u16, 14, 16] { + let top = ((1u32 << depth) - 1) as u16; + let mut px = image(w, h, 100, 4000); + px[5 * w + 9] = top; + let file = |px: &[u16]| format8_file(w, h, px, &[30, 34], &TABLE16, &|t| replace(t, BITS, Value::Short(vec![depth]))); + let r = crate::decode(&file(&px)).unwrap_or_else(|e| panic!("{depth}-bit: {e}")); + assert_eq!(samples(&r)[5 * w + 9], top, "{depth}-bit"); + if depth < 16 { + px[5 * w + 9] = top + 1; + let got = crate::decode(&file(&px)); + assert!(matches!(&got, Err(RawError::Corrupt(m)) if m.contains("outside")), "{depth}-bit: {got:?}"); + } + } + } + + #[test] + fn format8_tiles_in_two_rows_round_trip() { + // the DC-S1RM2 high-resolution mode: four columns by two rows of tiles, tag 0x0043 = 2 + let (w, h) = (96usize, 40usize); + let px = image(w, h, 520, 16000); + for (widths, rows) in [(vec![24, 24, 24, 24], vec![24, 16]), (vec![40, 56], vec![20, 20]), (vec![96], vec![8, 32])] { + let bytes = format8_tiled_file(w, &px, &widths, &rows, &TABLE14, &|_| {}); + let r = crate::decode(&bytes).unwrap_or_else(|e| panic!("{widths:?} x {rows:?}: {e}")); + assert_eq!(samples(&r), px, "{widths:?} x {rows:?}"); + assert_eq!(crate::probe_info(&bytes).unwrap(), r.info()); + } + // a tile row count that disagrees with the tiles, or tiles that leave a gap, are refused + let case = |tag: u16, v: Value| crate::decode(&format8_tiled_file(w, &px, &[48, 48], &[20, 20], &TABLE14, &|t| replace(t, tag, v.clone()))); + assert!(case(TILE_ROWS, Value::Short(vec![1])).is_err_and(|e| matches!(e, RawError::Unsupported(_)))); + assert!(case(STRIP_HEIGHTS, undefined(&[4, 20, 20, 20, 12])).is_err_and(|e| matches!(e, RawError::Corrupt(_)))); + } + + #[test] + fn format8_cell_is_coded_column_by_column() { + // One 2×2 cell coded by hand (not with `encode_strip`, which shares the decoder's site order): the values + // 1, 2, 3, 4 in stream order. Row order would put the two greens of an RGGB cell on each other's sites. + let mut bits: Vec = vec![]; + let mut put = |n: u32, v: u32| bits.extend((0..n).rev().map(|i| v >> i & 1 == 1)); + for v in [1u32, 2, 3, 4] { + let c = 32 - v.leading_zeros(); + let (len, code) = TABLE14[c as usize]; + put(len, code); + put(c, v); + } + let n = bits.len() as u64; + let mut data: Vec = bits.chunks(8).map(|b| b.iter().enumerate().fold(0u8, |a, (i, &x)| a | (x as u8) << i)).collect(); + data.extend([0; 8]); + let (mut top, mut bottom) = ([0u16; 2], [0u16; 2]); + decode_strip( + &Strip { data: &data, left: 0, top: 0, width: 2, height: 2, bits: n }, + &code_table(&TABLE14).unwrap(), + 16383, + vec![&mut top, &mut bottom], + ) + .unwrap(); + assert_eq!((top, bottom), ([1, 3], [2, 4]), "top-left, bottom-left, top-right, bottom-right"); + } + + fn replace(tags: &mut Vec<(u16, Value)>, tag: u16, v: Value) { + tags.retain(|(t, _)| *t != tag); + tags.push((tag, v)); + } + + fn undefined(v: &[u16]) -> Value { + Value::Undefined(v.iter().flat_map(|x| x.to_le_bytes()).collect()) + } + + #[test] + fn format8_rejects_what_it_cannot_decode() { + let (w, h) = (64usize, 30usize); + let px = image(w, h, 520, 16000); + let case = |tag: u16, v: Value| crate::decode(&format8_file(w, h, &px, &[32, 32], &TABLE14, &|t| replace(t, tag, v.clone()))); + // an auxiliary tag with a value no sample has + assert!(matches!(case(0x003b, Value::Short(vec![4095])), Err(RawError::Unsupported(_)))); + // strips that don't tile the width, or extend past the file + assert!(matches!(case(STRIP_WIDTHS, undefined(&[2, 32, 30])), Err(RawError::Corrupt(_)))); + assert!(matches!(case(STRIP_BITS, undefined(&[2, 0, 256, 0, 256])), Err(RawError::Corrupt(_)))); + // a strip whose stated length ends before its data does + assert!(matches!(case(STRIP_BITS, undefined(&[2, 64, 0, 64, 0])), Err(RawError::Corrupt(_)))); + // strips stating fewer bits than they have samples (one bit each at least): refused while reading the layout, + // so before the sensor-sized buffer is made (probing, which never decodes, refuses it too) + let short_bits = format8_file(w, h, &px, &[32, 32], &TABLE14, &|t| replace(t, STRIP_BITS, undefined(&[2, 959, 0, 959, 0]))); + assert!(matches!(crate::probe_info(&short_bits), Err(RawError::Corrupt(_)))); + assert!(matches!(crate::decode(&short_bits), Err(RawError::Corrupt(_)))); + // only the second strip short (the first keeps its real length), by one bit or entirely: still refused, naming it + let second_strip_bits = |stated: u32| { + format8_file(w, h, &px, &[32, 32], &TABLE14, &|t| { + let Some((_, Value::Undefined(b))) = t.iter_mut().find(|(tag, _)| *tag == STRIP_BITS) else { panic!("strip bits tag") }; + b[6..10].copy_from_slice(&[stated as u8, (stated >> 8) as u8, (stated >> 16) as u8, (stated >> 24) as u8]); + }) + }; + for stated in [959, 0] { + let file = second_strip_bits(stated); + for got in [crate::probe_info(&file).map(|_| ()), crate::decode(&file).map(|_| ())] { + assert!(matches!(&got, Err(RawError::Corrupt(m)) if m.contains("strip 1 states")), "{stated} bits: {got:?}"); + } + } + // code tables: two categories with the same code, a code longer than its length, missing entries + let table: Vec = std::iter::once(17).chain(TABLE14.iter().flat_map(|&(l, c)| [l as u16, c as u16])).collect(); + let mut dup = table.clone(); + dup[1..3].copy_from_slice(&[2, 0]); + assert!(matches!(case(HUFFMAN_CODES, undefined(&dup)), Err(RawError::Corrupt(_)))); + let mut long = table.clone(); + long[1..3].copy_from_slice(&[3, 9]); + assert!(matches!(case(HUFFMAN_CODES, undefined(&long)), Err(RawError::Corrupt(_)))); + assert!(matches!(case(HUFFMAN_CODES, undefined(&table[..20])), Err(RawError::Corrupt(_)))); + // the shadowed 12-bit codes of the 14-bit table are fine + assert!(case(HUFFMAN_CODES, undefined(&table)).is_ok()); + } + + #[test] + fn unsupported_and_truncated_files() { + let payload = [0u8; 64]; + for (format, bits) in [(1u16, 12u16), (3, 12), (4, 14), (7, 12), (9, 12)] { + assert!(matches!(crate::decode(&rw2(28, 4, bits, format, 128, &payload)), Err(RawError::Unsupported(_))), "{format}/{bits}"); + } + let tags = [short(BITS, 12), short(COMPRESSION, 34826)]; + let bytes = file(28, 4, &tags, &payload, &|off| vec![(t::STRIP_OFFSETS, Value::Long(vec![off]))]); + assert!(matches!(crate::decode(&bytes), Err(RawError::Unsupported(_)))); + for format in [4u16, 5, 6] { + let bytes = rw2(280, 20, 12, format, 128, &payload); + assert!(matches!(crate::decode(&bytes), Err(RawError::Corrupt(_))), "{format}"); + assert!(matches!(crate::probe_info(&bytes), Err(RawError::Corrupt(_))), "{format}"); + } + } +} diff --git a/crates/raw/src/shared.rs b/crates/raw/src/shared.rs index 65cb306a..63062f22 100644 --- a/crates/raw/src/shared.rs +++ b/crates/raw/src/shared.rs @@ -87,6 +87,119 @@ pub(crate) fn white_from_data(data: &[u16], bits: u32) -> f32 { } } +#[cfg(feature = "rw2")] +/// Output radii sampled to fit a radial warp. +const DISTORTION_SAMPLES: usize = 64; + +#[cfg(feature = "rw2")] +/// Least-squares radial warp from a camera's forward distortion map. +/// +/// `corrected` maps a source radius to the output radius it lands at; it must rise monotonically up to `t_max`, the +/// output radius of the image corners. Radii use the caller's unit; pass `WarpRectilinear`'s normalisation (distance +/// from the centre to the **farthest** corner = 1) so the result is directly a `WarpRectilinear` plane. The map is +/// scanned at 4096 points up to a source radius of `3 * t_max`: it must reach `t_max` within that range and rise at +/// every scan point (monotonicity is checked only there), otherwise `None`; a NaN at a scan point is rejected the +/// same way, other non-finite values by the fit. The map is inverted by bisection and +/// `source = u·(k0 + k1u² + k2u⁴ + k3u⁶)` is fitted over output radii `u·t_max`, `u` in [0, 1] ([`fit_radial_warp`]): +/// the returned `[k0, k1, k2, k3]` give the source radius in the caller's unit as a function of `u`. `None` too when +/// the worst fit error exceeds `max_fit_error`. +pub(crate) fn warp_from_forward_map(corrected: &dyn Fn(f64) -> f64, t_max: f64, max_fit_error: f64) -> Option<[f64; 4]> { + if !(t_max.is_finite() && t_max > 0.0) { + return None; + } + let mut top = None; + let mut prev = 0.0; + for i in 1..=4096 { + let r = 3.0 * t_max * i as f64 / 4096.0; + let v = corrected(r); + if v.is_nan() || v <= prev { + return None; + } + prev = v; + if v >= t_max { + top = Some(r); + break; + } + } + let top = top?; + let source = |t: f64| { + let (mut lo, mut hi) = (0.0, top); + for _ in 0..60 { + let m = 0.5 * (lo + hi); + if corrected(m) < t { lo = m } else { hi = m } + } + 0.5 * (lo + hi) + }; + fit_radial_warp(&|u| source(u * t_max), max_fit_error, 1.0) +} + +#[cfg(feature = "rw2")] +/// Least-squares `source = u·(k0 + k1u² + k2u⁴ + k3u⁶)` through `source_of_u` sampled at `u = i/64`, i = 0..=64 +/// (`u` = output radius, `WarpRectilinear` units, so the result is directly a plane's `[k0, k1, k2, k3]`). `None` when +/// a sample is not finite, the normal equations are singular, or the worst fit error at the samples up to the first +/// one at or beyond `gate_until` (the radii the image shows; 1 = all) exceeds `max_fit_error`. +pub(crate) fn fit_radial_warp(source_of_u: &dyn Fn(f64) -> f64, max_fit_error: f64, gate_until: f64) -> Option<[f64; 4]> { + let samples: Vec<(f64, f64)> = (0..=DISTORTION_SAMPLES) + .map(|i| { + let u = i as f64 / DISTORTION_SAMPLES as f64; + (u, source_of_u(u)) + }) + .collect(); + if samples.iter().any(|&(_, v)| !v.is_finite()) { + return None; + } + let basis = |u: f64| { + let u2 = u * u; + [u, u * u2, u * u2 * u2, u * u2 * u2 * u2] + }; + let mut a = [[0.0; 4]; 4]; + let mut rhs = [0.0; 4]; + for &(u, v) in &samples { + let f = basis(u); + for (i, fi) in f.iter().enumerate() { + for (j, fj) in f.iter().enumerate() { + a[i][j] += fi * fj; + } + rhs[i] += fi * v; + } + } + let k = solve4(a, rhs)?; + let fit = |u: f64| basis(u).iter().zip(&k).map(|(f, c)| f * c).sum::(); + // rounded up to a sample, so the band between the last sample inside and `gate_until` itself is checked too; + // capped at 1 (`min` also turns a NaN into 1: every sample is gated) + let n = DISTORTION_SAMPLES as f64; + let gate_until = (gate_until * n).ceil().min(n) / n; + let worst = samples.iter().filter(|&&(u, _)| u <= gate_until).map(|&(u, v)| (fit(u) - v).abs()).fold(0.0, f64::max); + (worst <= max_fit_error).then_some(k) +} + +#[cfg(feature = "rw2")] +/// Solve the 4 × 4 system `a·x = b` (Gaussian elimination with partial pivoting); `None` when singular or not finite. +fn solve4(mut a: [[f64; 4]; 4], mut b: [f64; 4]) -> Option<[f64; 4]> { + for c in 0..4 { + let p = (c..4).max_by(|&i, &j| a[i][c].abs().total_cmp(&a[j][c].abs()))?; + if a[p][c].is_nan() || a[p][c].abs() <= 1e-300 { + return None; + } + a.swap(c, p); + b.swap(c, p); + for r in c + 1..4 { + let f = a[r][c] / a[c][c]; + let (top, bottom) = a.split_at_mut(r); + for (x, p) in bottom[0][c..].iter_mut().zip(&top[c][c..]) { + *x -= f * p; + } + b[r] -= f * b[c]; + } + } + let mut x = [0.0; 4]; + for c in (0..4).rev() { + let s: f64 = (c + 1..4).map(|k| a[c][k] * x[k]).sum(); + x[c] = (b[c] - s) / a[c][c]; + } + x.iter().all(|v| v.is_finite()).then_some(x) +} + #[cfg(test)] pub(crate) fn declare_strip_len(file: &mut [u8], len: u32, declared: u32) { let big = file.starts_with(b"MM"); @@ -186,4 +299,55 @@ pub(crate) mod tests { let w = white_from_data(&d, 14); assert!((13500.0..=13583.0).contains(&w), "{w}"); } + + #[cfg(feature = "rw2")] + #[test] + fn warp_from_forward_map_recovers_identity_and_rejects_folds() { + let k = warp_from_forward_map(&|r| r, 1.0, 1e-9).unwrap(); + assert!((k[0] - 1.0).abs() < 1e-9 && k[1].abs() < 1e-9 && k[2].abs() < 1e-9 && k[3].abs() < 1e-9); + let barrel = |r: f64| r * (1.0 + 0.05 * r * r); // mild barrel: source < output radius + let t_max = 0.9; + let k = warp_from_forward_map(&barrel, t_max, 1e-3).unwrap(); + assert!(k[0] < t_max && k[0] > 0.8 && k[1] < 0.0); // source radius at u = 1 is below t_max + for i in 0..=20 { + let u = i as f64 / 20.0; + let s = u * (k[0] + u * u * (k[1] + u * u * (k[2] + u * u * k[3]))); + assert!((barrel(s) - u * t_max).abs() < 1e-5, "u {u}"); + } + assert!(warp_from_forward_map(&|r| if r < 0.5 { r } else { 0.5 - (r - 0.5) }, 1.0, 1e-3).is_none()); // fold, never reaches t_max + // a tiny dip (two scan points wide, 1e-4 deep) that still reaches t_max and fits well: only the monotonicity check rejects it + let dip = |r: f64| { + if r < 0.5 { + r + } else if r < 0.5015 { + 0.5 - 0.05 * (r - 0.5) + } else { + r - 0.001575 + } + }; + assert!(warp_from_forward_map(&dip, 1.0, 1e-2).is_none()); + assert!(warp_from_forward_map(&|_| f64::NAN, 1.0, 1e-3).is_none()); + } + + #[cfg(feature = "rw2")] + #[test] + fn fit_radial_warp_identity_and_gate() { + let k = fit_radial_warp(&|u| u, 1e-9, 1.0).unwrap(); + assert!((k[0] - 1.0).abs() < 1e-9 && k[1].abs() < 1e-9 && k[2].abs() < 1e-9 && k[3].abs() < 1e-9, "{k:?}"); + assert!(fit_radial_warp(&|u| if u > 0.5 { f64::NAN } else { u }, 1.0, 1.0).is_none(), "a NaN sample"); + assert!(fit_radial_warp(&|u| if u > 0.5 { f64::INFINITY } else { u }, 1.0, 1.0).is_none(), "an infinite sample"); + assert!(fit_radial_warp(&|u: f64| u.sqrt(), 1.0, 1.0).is_some()); + assert!(fit_radial_warp(&|u: f64| u.sqrt(), 0.0, 1.0).is_none(), "not a polynomial: the gate rejects it"); + // a step at u = 1 alone: worst fit error 0.057 over all samples, 0.006 over u <= 0.9 (the fit is the same) + let step = |u: f64| if u > 0.99 { 1.1 } else { u }; + assert!(fit_radial_warp(&step, 0.01, 1.0).is_none()); + assert_eq!(fit_radial_warp(&step, 0.01, 0.9), fit_radial_warp(&step, 1.0, 1.0), "gated only where shown"); + // the first sample at or beyond `gate_until` is gated too (the band up to it is shown): 0.023 up to + // u = 62/64, 0.032 at 63/64 + assert!(fit_radial_warp(&step, 0.03, 62.0 / 64.0).is_some(), "a sample exactly at the bound"); + assert!(fit_radial_warp(&step, 0.03, 0.98).is_none(), "63/64 is the next sample beyond 0.98"); + // capped at 1, and a NaN gates every sample (fails closed) + assert_eq!(fit_radial_warp(&step, 1.0, f64::INFINITY), fit_radial_warp(&step, 1.0, 1.0)); + assert!(fit_radial_warp(&step, 0.01, f64::NAN).is_none(), "a NaN gate"); + } } diff --git a/crates/raw/tests/it/vendor_robust.rs b/crates/raw/tests/it/vendor_robust.rs index 65fb139f..4df4912d 100644 --- a/crates/raw/tests/it/vendor_robust.rs +++ b/crates/raw/tests/it/vendor_robust.rs @@ -58,8 +58,119 @@ fn cr2() -> Vec { TiffWriter::new(ByteOrder::Little, false).write(&[ifd0, blank(), blank(), ifd3]).unwrap() } +/// A Panasonic RW2: sensor size, bit depth and raw format in IFD0, 12-bit packed blocks. +fn rw2() -> Vec { + let mut ifd = IfdBuilder::new(); + ifd.set(0x0002, Value::Short(vec![W as u16])); + ifd.set(0x0003, Value::Short(vec![H as u16])); + ifd.set(0x0004, Value::Short(vec![1])); + ifd.set(0x0005, Value::Short(vec![2])); + ifd.set(0x0006, Value::Short(vec![H as u16 - 1])); + ifd.set(0x0007, Value::Short(vec![W as u16 - 2])); + ifd.set(0x000a, Value::Short(vec![12])); + ifd.set(0x002d, Value::Short(vec![5])); + ifd.set(0x001c, Value::Short(vec![64])); + ifd.set(0x001d, Value::Short(vec![64])); + ifd.set(0x001e, Value::Short(vec![64])); + ifd.set(t::MAKE, Value::Ascii("Panasonic".into())); + let need = (W * H).div_ceil(10) as usize * 16; + ifd.set_image(ImageData::Strips { rows_per_strip: H, strips: vec![(0..need).map(|i| (i * 29) as u8).collect()] }); + let mut b = TiffWriter::new(ByteOrder::Little, false).write(&[ifd]).unwrap(); + b[..4].copy_from_slice(b"IIU\0"); + b +} + +/// A Panasonic RW2 of raw format `format` (or, with `format` 0, the oldest 16-bit-word kind): sensor size, borders, +/// bit depth, black levels, an in-camera crop and `n` bytes of payload at the end, at `RawDataOffset`. Format 8 also +/// gets two strips covering the payload and a code table. +fn rw2_of(format: u16, bits: u16, n: usize) -> Vec { + let u16s = |v: &[u16]| Value::Undefined(v.iter().flat_map(|x| x.to_le_bytes()).collect()); + let build = |off: u32| { + let mut ifd = IfdBuilder::new(); + for (tag, v) in [(2, W as u16), (3, H as u16), (4, 2), (5, 4), (6, H as u16), (7, W as u16), (9, 2), (0x0a, bits), (0x1c, 128), (0x1d, 128)] { + ifd.set(tag, Value::Short(vec![v])); + } + for (tag, v) in [(0x2f, 2u16), (0x30, 6), (0x31, H as u16 - 2), (0x32, W as u16 - 2)] { + ifd.set(tag, Value::Short(vec![v])); + } + if format == 0 { + ifd.set(0x000b, Value::Short(vec![34828])); + ifd.set(t::STRIP_OFFSETS, Value::Long(vec![off])); + } else { + ifd.set(0x002d, Value::Short(vec![format])); + ifd.set(0x0118, Value::Long(vec![off])); + } + if format == 8 { + let table = [ + (6, 62), + (7, 126), + (6, 61), + (5, 28), + (4, 12), + (3, 4), + (3, 2), + (2, 0), + (3, 3), + (3, 5), + (4, 13), + (5, 29), + (6, 60), + (8, 254), + (8, 255), + (12, 4094), + (12, 4095), + ]; + ifd.set(0x0040, u16s(&std::iter::once(17).chain(table.iter().flat_map(|&(l, c)| [l, c])).collect::>())); + let half = (n / 2) as u32; + let (a, b) = (off, off + half); + ifd.set(0x0044, u16s(&[2, a as u16, (a >> 16) as u16, b as u16, (b >> 16) as u16])); + ifd.set(0x0045, u16s(&[2, 0, 0, W as u16 / 2, 0])); + ifd.set(0x0046, u16s(&[2, (half * 8) as u16, ((half * 8) >> 16) as u16, (half * 8) as u16, ((half * 8) >> 16) as u16])); + ifd.set(0x0047, u16s(&[2, W as u16 / 2, W as u16 / 2])); + ifd.set(0x0048, u16s(&[2, H as u16, H as u16])); + } + ifd.set(t::MAKE, Value::Ascii("Panasonic".into())); + let mut b = TiffWriter::new(ByteOrder::Little, false).write(&[ifd]).unwrap(); + b[..4].copy_from_slice(b"IIU\0"); + b + }; + let off = build(0).len(); + let mut b = build(off as u32); + if format == 8 { + // a valid stream: the category-0 code `111110` over and over, each byte's bits from bit 0 up (every sample 0, + // inside any bit depth) + let bit = |k: usize| u8::from(k % 6 != 5); + b.extend((0..n).map(|i| (0..8).fold(0u8, |a, j| a | bit(i * 8 + j) << j))); + } else { + b.extend((0..n).map(|i| (i * 89 % 253) as u8)); + } + b +} + +/// One RW2 per encoding with a sample decoder of its own (formats 4, 6 at 14 and 12 bits, 8, and 16-bit words). +fn rw2_formats() -> Vec> { + vec![ + rw2_of(4, 12, (W * H).div_ceil(14) as usize * 16), + rw2_of(6, 14, (W * H).div_ceil(11) as usize * 16), + rw2_of(6, 12, (W * H).div_ceil(14) as usize * 16), + rw2_of(8, 14, (W * H) as usize * 2), + rw2_of(0, 12, (W * H) as usize * 2), + ] +} + +/// The RW2 samples decode as they are, so the mutations reach every format's sample decoder. +#[cfg(feature = "rw2")] +#[test] +fn rw2_samples_decode_unmutated() { + for s in &rw2_formats() { + let r = craft_raw::decode(s).unwrap(); + assert_eq!((r.width, r.height), (W as usize, H as usize)); + } +} + fn samples() -> Vec> { vec![ + rw2(), cr2(), cfa_tiff("NIKON CORPORATION", 1, 12, ByteOrder::Big), cfa_tiff("NIKON CORPORATION", 34713, 14, ByteOrder::Big), @@ -69,6 +180,9 @@ fn samples() -> Vec> { cfa_tiff("SAMSUNG", 32769, 14, ByteOrder::Big), cfa_tiff("ACME", 1, 16, ByteOrder::Little), ] + .into_iter() + .chain(rw2_formats()) + .collect() } fn exercise(bytes: &[u8]) {