diff --git a/crates/engine/src/camera_preview.rs b/crates/engine/src/camera_preview.rs index 54b308b13..3c8c74bae 100644 --- a/crates/engine/src/camera_preview.rs +++ b/crates/engine/src/camera_preview.rs @@ -31,7 +31,12 @@ pub(crate) fn fit_preview(raw: &RawImage, bytes: &[u8], transform: &CameraTransf } // Fixed, bounded proxy: the selected look cannot depend on thumbnail/export resolution. let k = (crop.width.max(crop.height).div_ceil(384).max(2)).div_ceil(2) * 2; - let sensor = raw.develop_binned(k, 0.99).ok()??; + let sensor = match raw.develop_binned(k, 0.99).ok()? { + Some(s) => s, + // reduced-size (M/S) ARWs are already RGB (issue #147): no demosaic, the box fit below does the binning + None if raw.cpp == 3 && raw.cfa.is_none() => raw.develop(lightcraft_raw::Method::Bilinear).ok()?, + None => return None, + }; let mut sensor = fit(&sensor, 96, 96, Filter::Box); let reference = fit(&reference, sensor.width, sensor.height, Filter::Box); let gain = 2f32.powf(transform.baseline_exposure as f32); diff --git a/crates/raw/examples/ycc_dump.rs b/crates/raw/examples/ycc_dump.rs new file mode 100644 index 000000000..60890a20c --- /dev/null +++ b/crates/raw/examples/ycc_dump.rs @@ -0,0 +1,46 @@ +//! TEMPORARY analysis helper (not committed): decode every tile of a subsampled lossless ARW into full planes. +//! Output: per plane `u32 w, u32 h` + `u16` samples, planes concatenated. +use lightcraft_tiff::Tiff; +use lightcraft_tiff::image::chunk_bytes; + +fn main() { + let a: Vec = std::env::args().skip(1).collect(); + let bytes = std::fs::read(&a[0]).unwrap(); + let tiff = Tiff::parse(&bytes).unwrap(); + let raw = tiff.all_ifds().into_iter().max_by_key(|i| i.u64(256).unwrap_or(0) * i.u64(257).unwrap_or(0)).unwrap(); + let info = raw.image().unwrap(); + let (w, h) = (info.width as usize, info.height as usize); + let chunks = info.chunks(bytes.len() as u64); + let mut planes: Vec<(usize, usize, Vec)> = Vec::new(); + for c in &chunks { + let src = chunk_bytes(&bytes, c).unwrap(); + let f = lightcraft_raw::ljpeg::decode_planar(src, 1 << 24).unwrap(); + if planes.is_empty() { + for p in &f.planes { + let (sx, sy) = (f.width / p.width, f.height / p.height); + planes.push((w / sx, h / sy, vec![0; (w / sx) * (h / sy)])); + } + } + for (k, p) in f.planes.iter().enumerate() { + let (sx, sy) = (f.width / p.width, f.height / p.height); + let (pw, ph, ref mut d) = planes[k]; + let (x0, y0) = (c.x as usize / sx, c.y as usize / sy); + for y in 0..p.height { + for x in 0..p.width { + if x0 + x < pw && y0 + y < ph { + d[(y0 + y) * pw + x0 + x] = p.data[y * p.width + x]; + } + } + } + } + } + let mut out = Vec::new(); + for (pw, ph, d) in &planes { + out.extend_from_slice(&(*pw as u32).to_le_bytes()); + out.extend_from_slice(&(*ph as u32).to_le_bytes()); + for v in d { + out.extend_from_slice(&v.to_le_bytes()); + } + } + std::fs::write(&a[1], out).unwrap(); +} diff --git a/crates/raw/src/ljpeg.rs b/crates/raw/src/ljpeg.rs index cb2d04675..98890c121 100644 --- a/crates/raw/src/ljpeg.rs +++ b/crates/raw/src/ljpeg.rs @@ -1,9 +1,11 @@ //! Lossless JPEG ("LJ92"): ITU-T T.81 / ISO 10918-1 process 14 (Annex H), Huffman-coded, as used by DNG, //! Canon CR2 and Sony lossless ARW. //! -//! - All seven predictors, 2–16 bit precision, 1–4 interleaved components (1×1 sampling), point transform, -//! restart intervals, `SSSS = 16` differences. -//! - [`decode`] returns the frame as interleaved samples in raster order (`width × height × components`), +//! - All seven predictors, 2–16 bit precision, 1–4 interleaved components, point transform, restart intervals, +//! `SSSS = 16` differences. +//! - Subsampled components (sampling factors other than 1×1, e.g. Sony's M/S-size lossless ARW with 4:2:0 / +//! 4:2:2 YCbCr): [`decode_planar`] returns one plane per component. +//! - [`decode`] returns a 1×1-sampled frame as interleaved samples in raster order (`width × height × components`), //! which callers re-tile (DNG tiles, CR2 slices). //! - [`encode`] writes a conforming stream with per-component optimal Huffman tables built with the //! T.81 Annex K.2 procedure — used by tests (bit-exact round trips) and later by DNG export. @@ -207,6 +209,8 @@ struct Header { height: usize, width: usize, comps: Vec, + /// Horizontal and vertical sampling factors per component (T.81 `Hi`, `Vi`). + sampling: Vec<(usize, usize)>, tables: [Option; 4], table_for: Vec, predictor: u8, @@ -228,6 +232,7 @@ fn parse_header(d: &[u8]) -> Result { height: 0, width: 0, comps: vec![], + sampling: vec![], tables: [None, None, None, None], table_for: vec![], predictor: 1, @@ -272,10 +277,16 @@ fn parse_header(d: &[u8]) -> Result { } for c in 0..n { let (id, samp) = (seg[6 + 3 * c], seg[7 + 3 * c]); - if samp != 0x11 { - return Err(RawError::Unsupported("lossless JPEG with subsampled components".into())); + let (hs, vs) = ((samp >> 4) as usize, (samp & 15) as usize); + if !(1..=4).contains(&hs) || !(1..=4).contains(&vs) { + return Err(err("bad sampling factor")); } h.comps.push(id); + h.sampling.push((hs, vs)); + } + // T.81 B.2.2: at most 10 data units per MCU + if h.sampling.iter().map(|(a, b)| a * b).sum::() > 10 { + return Err(err("too many samples per MCU")); } if !(2..=16).contains(&h.precision) { return Err(err("bad precision")); @@ -343,9 +354,18 @@ pub fn frame_info(d: &[u8]) -> Result<(usize, usize, usize, u8), RawError> { Ok((h.width, h.height, h.comps.len(), h.precision)) } -/// Decode a lossless JPEG stream. `max_samples` bounds the allocation (hostile headers). +/// The components' sampling factors `(H, V)` from the frame header (all `(1, 1)` unless subsampled). +pub fn sampling_factors(d: &[u8]) -> Result, RawError> { + Ok(parse_header(d)?.sampling) +} + +/// Decode a lossless JPEG stream. `max_samples` bounds the allocation (hostile headers). Subsampled frames +/// (sampling factors other than 1×1) are decoded by [`decode_planar`]. pub fn decode(d: &[u8], max_samples: usize) -> Result { let h = parse_header(d)?; + if h.sampling.iter().any(|&s| s != (1, 1)) { + return Err(RawError::Unsupported("lossless JPEG with subsampled components (planar decode only)".into())); + } let (w, ht, nc) = (h.width, h.height, h.comps.len()); if w == 0 || ht == 0 { return Err(err("zero frame size")); @@ -430,6 +450,163 @@ pub fn decode(d: &[u8], max_samples: usize) -> Result { Ok(Frame { width: w, height: ht, components: nc, precision: h.precision, predictor: pred, point_transform: h.pt, data: out }) } +/// One component of a [`PlanarFrame`]. +#[derive(Clone, Debug, PartialEq, Eq)] +pub struct Plane { + /// Sampling factors `(H, V)`. + pub sampling: (usize, usize), + /// Size in samples: `⌈X · H / Hmax⌉ × ⌈Y · V / Vmax⌉` (T.81 A.1.1). + pub width: usize, + pub height: usize, + /// `width × height` samples in raster order. + pub data: Vec, +} + +/// A decoded lossless-JPEG frame with one plane per component, for frames with subsampled components. +#[derive(Clone, Debug, PartialEq, Eq)] +pub struct PlanarFrame { + /// Frame size `X × Y` (the size of a component with the largest sampling factors). + pub width: usize, + pub height: usize, + pub precision: u8, + pub predictor: u8, + pub point_transform: u8, + pub planes: Vec, +} + +/// Decode a lossless JPEG stream into one plane per component, honouring sampling factors (ITU-T T.81 A.1.1, +/// A.2.3 and Annex H): the scan is a sequence of MCUs, each holding `Hi × Vi` samples of component `i` (in raster +/// order within the MCU) for every component in turn. Each sample is predicted from its neighbours in its own +/// component (H.1.2.1): the first sample of the scan or of a restart interval from `2^(P − Pt − 1)`, the rest of +/// that component line from the left neighbour, all samples not in the first column with the scan's predictor. +/// +/// First column, with vertical sampling `V > 1`: the first sample of each *MCU row* is predicted from the first +/// sample of the previous MCU row (`V` lines up), the first samples of the other lines of the MCU row from the +/// line above. T.81 reads as "the line above" for every line; the rule used here was established black-box on +/// Sony's 4:2:0 lossless ARWs (ILCE-7M4 M-size), where it is the only candidate that decodes without row-wide +/// steps (mean absolute step between adjacent row means 10.6 vs 45.6 for "the line above" on a 512-line tile, +/// and no wrapped samples, against 17 % wrapped over the whole image). With `V = 1` both readings agree. +/// +/// Restart intervals must span whole MCU rows. Works for 1×1 sampling too (then the planes are the +/// de-interleaved [`decode`] output). `max_samples` bounds the allocation. +pub fn decode_planar(d: &[u8], max_samples: usize) -> Result { + let h = parse_header(d)?; + let (w, ht, nc) = (h.width, h.height, h.comps.len()); + if w == 0 || ht == 0 { + return Err(err("zero frame size")); + } + let hmax = h.sampling.iter().map(|s| s.0).max().unwrap_or(1); + let vmax = h.sampling.iter().map(|s| s.1).max().unwrap_or(1); + let (mcux, mcuy) = (w.div_ceil(hmax), ht.div_ceil(vmax)); + // planes are decoded at whole-MCU size and trimmed afterwards + let padded: Vec<(usize, usize)> = h.sampling.iter().map(|&(hs, vs)| (mcux * hs, mcuy * vs)).collect(); + let total = padded + .iter() + .try_fold(0usize, |acc, &(pw, ph)| pw.checked_mul(ph).and_then(|n| acc.checked_add(n))) + .ok_or_else(|| err("frame too large"))?; + if total > max_samples { + return Err(RawError::Limit("lossless JPEG frame larger than expected")); + } + let entropy = &d[h.scan_start..]; + if (entropy.len() as u64 + 64) * 8 < total as u64 { + return Err(err("entropy data too short for frame")); + } + if h.restart > 0 && h.restart % mcux != 0 { + return Err(RawError::Unsupported("lossless JPEG restart interval inside an MCU row of a subsampled frame".into())); + } + let tables: Vec<&Huffman> = h + .table_for + .iter() + .map(|&t| h.tables.get(t).and_then(Option::as_ref).ok_or_else(|| err("scan uses an undefined Huffman table"))) + .collect::>()?; + let mut planes: Vec> = padded.iter().map(|&(pw, ph)| vec![0u16; pw * ph]).collect(); + let mut br = BitReader::new(entropy); + let init = 1i32 << (h.precision - h.pt - 1); + let pred = h.predictor; + let mut restart_left = h.restart; + // MCU row where the scan / current restart interval started: its first line per component is the "first line" + let mut interval_row = 0usize; + let mut reset_pending = vec![true; nc]; + for my in 0..mcuy { + for mx in 0..mcux { + if h.restart > 0 { + if restart_left == 0 { + if !br.restart() { + return Err(err("missing restart marker")); + } + restart_left = h.restart; + interval_row = my; + reset_pending.iter_mut().for_each(|r| *r = true); + } + restart_left -= 1; + } + for c in 0..nc { + let (hs, vs) = h.sampling[c]; + let pw = padded[c].0; + let plane = &mut planes[c]; + let first_line = interval_row * vs; + for by in 0..vs { + let y = my * vs + by; + for bx in 0..hs { + let x = mx * hs + bx; + let i = y * pw + x; + let p = if reset_pending[c] { + reset_pending[c] = false; + init + } else if y == first_line { + // x > 0 here: the line's first sample was the reset sample + plane[i - 1] as i32 + } else if pred == 0 { + 0 + } else if x == 0 { + // the first sample of an MCU row: from the first sample of the previous MCU row + // (see the function docs); the first sample of any other line: from the one above + if by == 0 { plane[i - vs * pw] as i32 } else { plane[i - pw] as i32 } + } else { + let ra = plane[i - 1] as i32; + let rb = plane[i - pw] as i32; + let rc = plane[i - pw - 1] as i32; + match pred { + 1 => ra, + 2 => rb, + 3 => rc, + 4 => ra + rb - rc, + 5 => ra + ((rb - rc) >> 1), + 6 => rb + ((ra - rc) >> 1), + _ => (ra + rb) >> 1, + } + }; + let ssss = tables[c].decode(&mut br)?; + let diff = diff_value(&mut br, ssss); + plane[i] = ((p + diff) & 0xffff) as u16; + } + } + } + } + if br.overrun > 64 * 8 { + return Err(err("entropy data exhausted")); + } + } + let planes = planes + .into_iter() + .zip(&h.sampling) + .zip(&padded) + .map(|((mut data, &(hs, vs)), &(pw, _))| { + let (cw, ch) = ((w * hs).div_ceil(hmax), (ht * vs).div_ceil(vmax)); + if cw != pw { + data = data.chunks_exact(pw).take(ch).flat_map(|row| row[..cw].iter().copied()).collect(); + } else { + data.truncate(cw * ch); + } + if h.pt > 0 { + data.iter_mut().for_each(|v| *v <<= h.pt); + } + Plane { sampling: (hs, vs), width: cw, height: ch, data } + }) + .collect(); + Ok(PlanarFrame { width: w, height: ht, precision: h.precision, predictor: pred, point_transform: h.pt, planes }) +} + // ------------------------------------------------------------------------------------------------ encoder /// Build Huffman code lengths from frequencies (T.81 Annex K.2, Figures K.1–K.3), limited to 16 bits. @@ -700,6 +877,143 @@ pub fn encode(data: &[u16], width: usize, height: usize, components: usize, prec out } +/// Test encoder for subsampled frames: `planes[c]` is component `c` at `⌈w · H / Hmax⌉ × ⌈h · V / Vmax⌉` samples, +/// coded in interleaved MCUs with the prediction rules of [`decode_planar`] (one Huffman table per component). +#[cfg(test)] +pub(crate) fn encode_planar( + planes: &[Vec], + w: usize, + h: usize, + sampling: &[(usize, usize)], + precision: u8, + predictor: u8, + restart: usize, +) -> Vec { + let nc = planes.len(); + let (hmax, vmax) = (sampling.iter().map(|s| s.0).max().unwrap_or(1), sampling.iter().map(|s| s.1).max().unwrap_or(1)); + let (mcux, mcuy) = (w.div_ceil(hmax), h.div_ceil(vmax)); + // pad every plane to whole MCUs by repeating its last column / row + let padded: Vec<(usize, Vec)> = planes + .iter() + .zip(sampling) + .map(|(p, &(hs, vs))| { + let (cw, ch) = ((w * hs).div_ceil(hmax), (h * vs).div_ceil(vmax)); + let (pw, ph) = (mcux * hs, mcuy * vs); + (pw, (0..pw * ph).map(|i| p[(i / pw).min(ch - 1) * cw + (i % pw).min(cw - 1)]).collect()) + }) + .collect(); + let init = 1i32 << (precision - 1); + // differences in coding order, with the component of each + let mut diffs: Vec<(usize, i32)> = Vec::new(); + let mut marks = Vec::new(); + let mut restart_left = restart; + let mut interval_row = 0; + let mut reset = vec![true; nc]; + for my in 0..mcuy { + for mx in 0..mcux { + if restart > 0 { + if restart_left == 0 { + marks.push(diffs.len()); + restart_left = restart; + interval_row = my; + reset.iter_mut().for_each(|r| *r = true); + } + restart_left -= 1; + } + for c in 0..nc { + let (hs, vs) = sampling[c]; + let (pw, ref p) = padded[c]; + for by in 0..vs { + for bx in 0..hs { + let (x, y) = (mx * hs + bx, my * vs + by); + let i = y * pw + x; + let pred = if reset[c] { + reset[c] = false; + init + } else if y == interval_row * vs { + p[i - 1] as i32 + } else if x == 0 { + if by == 0 { p[i - vs * pw] as i32 } else { p[i - pw] as i32 } + } else { + let (ra, rb, rc) = (p[i - 1] as i32, p[i - pw] as i32, p[i - pw - 1] as i32); + match predictor { + 1 => ra, + 2 => rb, + 3 => rc, + 4 => ra + rb - rc, + 5 => ra + ((rb - rc) >> 1), + 6 => rb + ((ra - rc) >> 1), + _ => (ra + rb) >> 1, + } + }; + let mut d = (p[i] as i32 - pred) & 0xffff; + if d > 32768 { + d -= 65536; + } + diffs.push((c, d)); + } + } + } + } + } + let tables: Vec<([u8; 16], Vec)> = (0..nc) + .map(|c| { + let mut freq = [0u32; 17]; + diffs.iter().filter(|(k, _)| *k == c).for_each(|(_, d)| freq[ssss_of(*d) as usize] += 1); + code_lengths(&freq) + }) + .collect(); + let mut out = vec![0xff, 0xd8, 0xff, 0xc3]; + out.extend_from_slice(&((8 + 3 * nc) as u16).to_be_bytes()); + out.push(precision); + out.extend_from_slice(&(h as u16).to_be_bytes()); + out.extend_from_slice(&(w as u16).to_be_bytes()); + out.push(nc as u8); + for (c, &(hs, vs)) in sampling.iter().enumerate() { + out.extend_from_slice(&[c as u8 + 1, (hs << 4 | vs) as u8, 0]); + } + for (c, (counts, values)) in tables.iter().enumerate() { + out.extend_from_slice(&[0xff, 0xc4]); + out.extend_from_slice(&((2 + 17 + values.len()) as u16).to_be_bytes()); + out.push(c as u8); + out.extend_from_slice(counts); + out.extend_from_slice(values); + } + if restart > 0 { + out.extend_from_slice(&[0xff, 0xdd, 0, 4]); + out.extend_from_slice(&(restart as u16).to_be_bytes()); + } + out.extend_from_slice(&[0xff, 0xda]); + out.extend_from_slice(&((6 + 2 * nc) as u16).to_be_bytes()); + out.push(nc as u8); + for c in 0..nc { + out.extend_from_slice(&[c as u8 + 1, (c as u8) << 4]); + } + out.extend_from_slice(&[predictor, 0, 0]); + let codes: Vec<[(u32, u8); 17]> = tables.iter().map(|(c, v)| canonical_codes(c, v)).collect(); + let mut bw = BitWriter { out, acc: 0, n: 0 }; + let mut rst = 0u8; + for (i, &(c, d)) in diffs.iter().enumerate() { + if marks.contains(&i) { + bw.flush(); + bw.out.extend_from_slice(&[0xff, 0xd0 + rst]); + rst = (rst + 1) & 7; + bw.acc = 0; + } + let s = ssss_of(d); + let (code, len) = codes[c][s as usize]; + bw.put(code, len as u32); + if s > 0 && s < 16 { + let v = if d < 0 { d - 1 } else { d }; + bw.put(v as u32 & ((1 << s) - 1), s as u32); + } + } + bw.flush(); + let mut out = bw.out; + out.extend_from_slice(&[0xff, 0xd9]); + out +} + #[cfg(test)] mod tests { use super::*; @@ -771,6 +1085,92 @@ mod tests { let _ = decode(&bad, usize::MAX); } + /// Smooth test planes for a `w × h` frame with the given sampling factors. + fn planar_input(w: usize, h: usize, sampling: &[(usize, usize)], seed: u64) -> Vec> { + let (hmax, vmax) = (sampling.iter().map(|s| s.0).max().unwrap(), sampling.iter().map(|s| s.1).max().unwrap()); + sampling + .iter() + .enumerate() + .map(|(c, &(hs, vs))| { + let (pw, ph) = ((w * hs).div_ceil(hmax), (h * vs).div_ceil(vmax)); + let n = noise(pw * ph, 12, seed + c as u64); + (0..pw * ph).map(|i| 1000 + (i / pw * 7 + i % pw * 3) as u16 % 2000 + n[i] % 300).collect() + }) + .collect() + } + + #[test] + fn planar_roundtrip_subsampled() { + for (w, h, sampling) in [ + (16usize, 12usize, vec![(2usize, 2usize), (1, 1), (1, 1)]), // 4:2:0 (Sony M-size) + (16, 12, vec![(2, 1), (1, 1), (1, 1)]), // 4:2:2 (Sony S-size) + (15, 11, vec![(2, 2), (1, 1), (1, 1)]), // partial MCUs at the right and bottom + (9, 7, vec![(1, 1), (1, 1)]), // no subsampling + (8, 8, vec![(1, 2), (2, 1), (1, 1)]), + ] { + let planes = planar_input(w, h, &sampling, 3); + for pred in 1..=7 { + let enc = encode_planar(&planes, w, h, &sampling, 15, pred, 0); + let f = decode_planar(&enc, usize::MAX).unwrap_or_else(|e| panic!("{w}×{h} {sampling:?} p{pred}: {e}")); + assert_eq!((f.width, f.height, f.predictor), (w, h, pred)); + for (k, p) in f.planes.iter().enumerate() { + assert_eq!(p.sampling, sampling[k]); + assert_eq!(p.data, planes[k], "{w}×{h} {sampling:?} predictor {pred} component {k}"); + } + } + } + // restart intervals of whole MCU rows + let s = [(2, 2), (1, 1), (1, 1)]; + let planes = planar_input(16, 12, &s, 9); + let enc = encode_planar(&planes, 16, 12, &s, 15, 1, 16); + assert_eq!(decode_planar(&enc, usize::MAX).unwrap().planes.into_iter().map(|p| p.data).collect::>(), planes); + // ... others are reported, not mis-decoded + let enc = encode_planar(&planes, 16, 12, &s, 15, 1, 5); + assert!(matches!(decode_planar(&enc, usize::MAX), Err(RawError::Unsupported(_)))); + // the interleaved decoder hands subsampled frames to the planar one + let enc = encode_planar(&planes, 16, 12, &s, 15, 1, 0); + assert!(matches!(decode(&enc, usize::MAX), Err(RawError::Unsupported(_)))); + assert_eq!(sampling_factors(&enc).unwrap(), s.to_vec()); + } + + #[test] + fn planar_first_column_follows_the_mcu_row() { + // a frame whose first column steps between MCU rows: decoding with "the line above" for the first sample + // of an MCU row would shift every later row; our encoder mirrors the rule of `decode_planar` + let s = [(2, 2), (1, 1), (1, 1)]; + let mut planes = planar_input(8, 8, &s, 5); + for (y, row) in planes[0].chunks_mut(8).enumerate() { + row[0] = 4000 + 500 * y as u16; + } + let enc = encode_planar(&planes, 8, 8, &s, 15, 1, 0); + assert_eq!(decode_planar(&enc, usize::MAX).unwrap().planes[0].data, planes[0]); + } + + #[test] + fn planar_rejects_bad_streams() { + let s = [(2, 2), (1, 1), (1, 1)]; + let planes = planar_input(16, 12, &s, 1); + let enc = encode_planar(&planes, 16, 12, &s, 15, 1, 0); + assert!(decode_planar(&[], 100).is_err()); + assert!(matches!(decode_planar(&enc, 100), Err(RawError::Limit(_)))); + for n in 0..enc.len() { + let _ = decode_planar(&enc[..n], usize::MAX); + } + let mut bad = enc.clone(); + for i in (20..bad.len()).step_by(3) { + bad[i] ^= 0x5a; + } + let _ = decode_planar(&bad, usize::MAX); + // sampling factors 0 or > 4, and more than 10 samples per MCU + let sof = enc.windows(2).position(|w| w == [0xff, 0xc3]).unwrap(); + for samp in [0x00u8, 0x50, 0x05, 0x44] { + let mut b = enc.clone(); + b[sof + 11] = samp; // first component's sampling byte + b[sof + 14] = samp; + assert!(decode_planar(&b, usize::MAX).is_err(), "sampling {samp:#x}"); + } + } + #[test] fn code_lengths_are_limited_to_16() { // Fibonacci-like frequencies produce very long codes before adjustment diff --git a/crates/raw/src/vendor/arw.rs b/crates/raw/src/vendor/arw.rs index fc019743e..34a443fa9 100644 --- a/crates/raw/src/vendor/arw.rs +++ b/crates/raw/src/vendor/arw.rs @@ -21,8 +21,9 @@ //! default 512 on 1″-sensor bodies such as the RX100 series), else 512 (14-bit) / 128 (12-bit). //! //! Also: uncompressed 16-bit ARW, and lossless-compressed ARW (Compression 7, ILCE-7M4 and later): LJ92 tiles whose -//! frames hold one 2×2 CFA cell per four-component sample ([`read_quad_tiles`]). Other lossless-JPEG layouts go -//! through the generic TIFF path. +//! frames hold one 2×2 CFA cell per four-component sample ([`read_quad_tiles`]); reduced-size (M/S) lossless ARWs +//! store subsampled YCbCr instead ([`decode_ycbcr`], issue #147). Other lossless-JPEG layouts go through the +//! generic TIFF path. use crate::tiffraw::{Packing, check_image, read_image_in}; use crate::{BlackLevel, Cfa, ColorData, Mode, OpcodeLists, RawData, RawError, RawFormat, RawImage, Rect, Result, ljpeg}; @@ -300,6 +301,199 @@ fn read_quad_tiles(bytes: &[u8], info: &ImageInfo) -> Result> { Ok(out) } +/// Black offset of the luma in Sony's YCbCr ARWs, in the samples' own units. The files' `ReferenceBlackWhite` +/// says 0, but on the CC0 samples the darkest pixels sit at ≈1060 and, registered against the same scene shot as +/// a full-size (CFA) ARW (ILCE-7RM5), luma = 0.997 × (white-balanced, black-subtracted CFA luma) + 1034. +const YCC_BLACK: f32 = 1024.0; +/// Sony's clip level of the white-balanced RGB behind the YCbCr samples: clipped highlights are stored as a +/// neutral luma plateau at 18 876 (ILCE-7M4 M-size sample), i.e. 17 852 above [`YCC_BLACK`]. +const YCC_WHITE: f32 = 18876.0; +/// The values written for the camera RGB we reconstruct: `YCC_OUT_BLACK` + black-subtracted camera RGB (a raw-like +/// scale), white level `YCC_OUT_BLACK + YCC_WHITE − YCC_BLACK` for every channel. +const YCC_OUT_BLACK: f32 = 512.0; +/// TIFF 6.0 `YCbCrCoefficients` (`LumaRed`, `LumaGreen`, `LumaBlue`). +const YCBCR_COEFFICIENTS: u16 = 0x0211; + +/// Sony's reduced-size lossless ARWs ("M" and "S" RAW size on the ILCE-7M4, -7RM5, -1 and later): TIFF +/// compression 7 with `PhotometricInterpretation` 6 (YCbCr), three samples per pixel, `YCbCrSubSampling` 2×2 (M) +/// or 2×1 (S). Each tile is one lossless-JPEG frame of the tile's size whose luma has sampling factors 2×2 / 2×1 +/// and whose two chroma components are 1×1 ([`ljpeg::decode_planar`]). The image is already demosaiced and white +/// balanced: luma and chroma follow the TIFF 6.0 YCbCr definition (section 21) with the file's +/// `YCbCrCoefficients` (BT.601: 0.299, 0.587, 0.114), chroma centred on 16 384 and on the luma's scale — checked +/// by registering an M- and an S-size ILCE-7RM5 file against the L-size (CFA) one of the same scene: luma +/// = 0.997·Y + 1034 (M) and chroma = 1.01 × the TIFF formula (S: Cb 1.014, Cr 1.008; M: Cb 1.006, Cr 0.89, the +/// least reliable of the four with the M frame registered at a non-integer scale). We invert that (with chroma +/// co-sited, `YCbCrPositioning` 2, and interpolated bilinearly), remove [`YCC_BLACK`] and divide by the file's white balance (`0x7313`), so +/// the result is linear camera RGB like a `LinearRaw` DNG, with the as-shot multipliers reported as usual. +fn decode_ycbcr( + bytes: &[u8], + mode: Mode, + tiff: &Tiff, + raw: &Ifd, + info: &ImageInfo, + wb: Option<[f32; 3]>, + mn: Option<&makernote::MakerNote>, +) -> Result { + let (w, h) = (info.width as usize, info.height as usize); + if info.samples_per_pixel != 3 || !matches!(info.layout, Layout::Tiles { .. }) { + return Err(RawError::Unsupported("Sony YCbCr ARW that is not tiled with three samples per pixel".into())); + } + check_image(bytes, info)?; + let coeffs = match raw.f64s(YCBCR_COEFFICIENTS).as_deref() { + Some([r, g, b]) if [r, g, b].iter().all(|v| (0.01..1.0).contains(*v)) && (r + g + b - 1.0).abs() < 0.01 => [*r as f32, *g as f32, *b as f32], + _ => [0.299, 0.587, 0.114], + }; + let gains = wb.unwrap_or([1.0; 3]); + let data = match mode { + Mode::Header => Vec::new(), + Mode::Full => { + let (luma, chroma, (hs, vs)) = read_ycbcr_planes(bytes, info)?; + ycbcr_to_camera_rgb(&luma, &chroma, w, h, (hs, vs), coeffs, gains) + } + }; + // one white level for all channels, as for a mosaic: Sony clipped after white balance, so a clipped highlight + // (R, G, B = clip / gains in camera RGB) becomes exactly white once the pipeline applies the multipliers + let white = vec![YCC_OUT_BLACK + (YCC_WHITE - YCC_BLACK)]; + let crop = match (raw.u64s(CROP_TOP_LEFT).as_deref(), raw.u64s(CROP_SIZE).as_deref()) { + (Some([x, y]), Some([cw, ch])) if *cw > 0 && *ch > 0 => Rect::new(*x as usize, *y as usize, *cw as usize, *ch as usize).clipped(w, h), + _ => default_crop(raw, mn, w, h), + }; + let mut metadata = lightcraft_meta::from_tiff(tiff); + metadata.width = Some(crop.width as u32); + metadata.height = Some(crop.height as u32); + let img = RawImage { + format: RawFormat::Arw, + width: w, + height: h, + cpp: 3, + data: RawData::U16(data), + cfa: None, + bits: info.bits() as u32, + black: BlackLevel::uniform(YCC_OUT_BLACK), + white, + active_area: Rect::new(0, 0, w, h), + crop, + orientation: Orientation::from_exif(tiff.ifds[0].u16(t::ORIENTATION).unwrap_or(1)), + color: ColorData::default(), + wb_multipliers: wb, + linearized: false, + opcodes: OpcodeLists::default(), + metadata, + }; + img.validate_for(mode)?; + Ok(img) +} + +/// Decode every tile of a subsampled YCbCr ARW into a full-size luma plane and two chroma planes (Cb, Cr) of +/// `⌈w / H⌉ × ⌈h / V⌉` samples; returns them with the luma sampling factors `(H, V)`. +type YccPlanes = (Vec, [Vec; 2], (usize, usize)); +fn read_ycbcr_planes(bytes: &[u8], info: &ImageInfo) -> Result { + let (w, h) = (info.width as usize, info.height as usize); + let chunks = info.chunks(bytes.len() as u64); + let decoded: Vec> = chunks + .par_iter() + .map(|c| { + let src = chunk_bytes(bytes, c).ok_or_else(|| RawError::Corrupt("tile offset past end of file".into()))?; + let f = ljpeg::decode_planar(src, (c.width as usize * c.height as usize * 3).max(1 << 16))?; + let ok = f.planes.len() == 3 + && (f.width, f.height) == (c.width as usize, c.height as usize) + && f.planes[1..].iter().all(|p| p.sampling == (1, 1)) + && matches!(f.planes[0].sampling, (1 | 2, 1 | 2)); + if !ok { + return Err(RawError::Unsupported(format!( + "Sony YCbCr ARW tile: {}×{} frame with sampling {:?}", + f.width, + f.height, + f.planes.iter().map(|p| p.sampling).collect::>() + ))); + } + Ok((*c, f)) + }) + .collect(); + let mut sampling = None; + let mut luma = vec![0u16; w * h]; + let mut chroma: [Vec; 2] = [Vec::new(), Vec::new()]; + let (mut cw, mut ch) = (0, 0); + let mut ok = 0usize; + let mut first_err = None; + for r in decoded { + let (c, f) = match r { + Ok(v) => v, + Err(e) => { + first_err.get_or_insert(e); + continue; + } + }; + let (hs, vs) = f.planes[0].sampling; + match sampling { + None => { + sampling = Some((hs, vs)); + (cw, ch) = (w.div_ceil(hs), h.div_ceil(vs)); + chroma = [vec![16384u16; cw * ch], vec![16384u16; cw * ch]]; + } + Some(s) if s != (hs, vs) => { + first_err.get_or_insert(RawError::Corrupt("Sony YCbCr ARW: tiles with different sampling".into())); + continue; + } + Some(_) => {} + } + ok += 1; + let (x0, y0) = (c.x as usize, c.y as usize); + let p = &f.planes[0]; + for y in 0..p.height.min(h.saturating_sub(y0)) { + let n = p.width.min(w.saturating_sub(x0)); + luma[(y0 + y) * w + x0..][..n].copy_from_slice(&p.data[y * p.width..][..n]); + } + let (cx0, cy0) = (x0 / hs, y0 / vs); + for (k, plane) in f.planes[1..].iter().enumerate() { + for y in 0..plane.height.min(ch.saturating_sub(cy0)) { + let n = plane.width.min(cw.saturating_sub(cx0)); + chroma[k][(cy0 + y) * cw + cx0..][..n].copy_from_slice(&plane.data[y * plane.width..][..n]); + } + } + } + match sampling { + Some(s) if ok > 0 => Ok((luma, chroma, s)), + _ => Err(first_err.unwrap_or_else(|| RawError::Corrupt("no decodable tiles".into()))), + } +} + +/// YCbCr → camera RGB (see [`decode_ycbcr`]): interleaved `w × h × 3` samples on a raw-like scale. +fn ycbcr_to_camera_rgb(luma: &[u16], chroma: &[Vec; 2], w: usize, h: usize, (hs, vs): (usize, usize), coeffs: [f32; 3], gains: [f32; 3]) -> Vec { + let (cw, ch) = (w.div_ceil(hs), h.div_ceil(vs)); + let [kr, kg, kb] = coeffs; + let inv = gains.map(|g| if g > 0.0 && g.is_finite() { 1.0 / g } else { 1.0 }); + let mut out = vec![0u16; w * h * 3]; + out.par_chunks_mut(w * 3).enumerate().for_each(|(y, row)| { + // co-sited chroma: sample (i, j) sits on luma (i·H, j·V) + let fy = y as f32 / vs as f32; + let j0 = (fy as usize).min(ch - 1); + let j1 = (j0 + 1).min(ch - 1); + let ty = fy - j0 as f32; + for x in 0..w { + let fx = x as f32 / hs as f32; + let i0 = (fx as usize).min(cw - 1); + let i1 = (i0 + 1).min(cw - 1); + let tx = fx - i0 as f32; + let at = |p: &Vec| { + let v = |i: usize, j: usize| f32::from(p.get(j * cw + i).copied().unwrap_or(16384)); + let top = v(i0, j0) + (v(i1, j0) - v(i0, j0)) * tx; + let bottom = v(i0, j1) + (v(i1, j1) - v(i0, j1)) * tx; + top + (bottom - top) * ty - 16384.0 + }; + let (cb, cr) = (at(&chroma[0]), at(&chroma[1])); + let yv = f32::from(luma.get(y * w + x).copied().unwrap_or(0)) - YCC_BLACK; + let r = cr * (2.0 - 2.0 * kr) + yv; + let b = cb * (2.0 - 2.0 * kb) + yv; + let g = (yv - kb * b - kr * r) / kg; + for (c, v) in [r, g, b].into_iter().enumerate() { + row[x * 3 + c] = (YCC_OUT_BLACK + v * inv[c]).round().clamp(0.0, 65535.0) as u16; + } + } + }); + out +} + /// `SR2Private` tags (ExifTool "Sony SR2Private" table; the IFD is referenced by IFD0's `DNGPrivateData`). const SR2_SUBIFD_OFFSET: u16 = 0x7200; const SR2_SUBIFD_LENGTH: u16 = 0x7201; @@ -383,7 +577,23 @@ pub(crate) fn decode(bytes: &[u8], mode: Mode) -> Result { if w * h > crate::MAX_SAMPLES { return Err(RawError::Limit("image too large")); } + let model = ifd0.string(t::MODEL).unwrap_or_default(); + let mn = tiff + .exif() + .and_then(|e| e.get(t::MAKER_NOTE)) + .and_then(|e| makernote::parse_makernote(bytes, e.offset, e.count() as u64, tiff.order, &ifd0.string(t::MAKE).unwrap_or_default())); + let wb = raw + .f64s(WB_RGGB) + .filter(|v| v.len() == 4 && v[1] > 0.0 && v[0] > 0.0 && v[3] > 0.0) + .map(|v| { + let g = (v[1] + v[2]) / 2.0; + [(v[0] / g) as f32, 1.0, (v[3] / g) as f32] + }) + .or_else(|| mn.as_ref().and_then(|m| tag2010_wb(&model, m.ifd.bytes(MN_TAG2010)?, m.order))); let bits = info.bits() as u32; + if info.compression == 7 && info.photometric == photometric::YCBCR { + return decode_ycbcr(bytes, mode, &tiff, raw, &info, wb, mn.as_ref()); + } let chunks = info.chunks(bytes.len() as u64); let strip_len: u64 = chunks.iter().map(|c| c.len).sum(); let (data, out_bits) = match info.compression { @@ -440,19 +650,6 @@ pub(crate) fn decode(bytes: &[u8], mode: Mode) -> Result { _ => BlackLevel::uniform(sr2_black(bytes, ifd0, tiff.order).filter(|_| scale_bits >= 14).unwrap_or(default_black)), }; let white = raw.f64(t::WHITE_LEVEL).map(|v| v as f32).filter(|v| *v > 0.0).unwrap_or_else(|| super::white_from_data(samples, scale_bits)); - let model = ifd0.string(t::MODEL).unwrap_or_default(); - let mn = tiff - .exif() - .and_then(|e| e.get(t::MAKER_NOTE)) - .and_then(|e| makernote::parse_makernote(bytes, e.offset, e.count() as u64, tiff.order, &ifd0.string(t::MAKE).unwrap_or_default())); - let wb = raw - .f64s(WB_RGGB) - .filter(|v| v.len() == 4 && v[1] > 0.0 && v[0] > 0.0 && v[3] > 0.0) - .map(|v| { - let g = (v[1] + v[2]) / 2.0; - [(v[0] / g) as f32, 1.0, (v[3] / g) as f32] - }) - .or_else(|| mn.as_ref().and_then(|m| tag2010_wb(&model, m.ifd.bytes(MN_TAG2010)?, m.order))); let crop = match (raw.u64s(CROP_TOP_LEFT).as_deref(), raw.u64s(CROP_SIZE).as_deref()) { (Some([x, y]), Some([cw, ch])) if *cw > 0 && *ch > 0 => Rect::new(*x as usize, *y as usize, *cw as usize, *ch as usize).clipped(w, h), _ => default_crop(raw, mn.as_ref(), w, h), @@ -602,6 +799,81 @@ mod tests { assert_eq!(sr2_black_levels(&[], Little), None); } + /// Camera RGB → Sony's YCbCr samples (the model `decode_ycbcr` inverts). + fn to_ycc(rgb: [f32; 3], gains: [f32; 3]) -> [f32; 3] { + let [r, g, b] = [rgb[0] * gains[0], rgb[1] * gains[1], rgb[2] * gains[2]]; + let y = 0.299 * r + 0.587 * g + 0.114 * b; + [y + YCC_BLACK, (b - y) / 1.772 + 16384.0, (r - y) / 1.402 + 16384.0] + } + + #[test] + fn ycbcr_tiles_decode_to_camera_rgb() { + // 4:2:0 (M) and 4:2:2 (S) images of 20×12 in 8×8 tiles, from camera RGB that is linear in x and y (so the + // bilinear interpolation of co-sited chroma is exact) and a white balance of 2.4 / 1.7 + let gains = [2.4f32, 1.0, 1.7]; + for (hs, vs) in [(2usize, 2usize), (2, 1)] { + let (w, h, tw) = (20usize, 12usize, 8usize); + let rgb_at = |x: usize, y: usize| [300.0 + 20.0 * x as f32, 900.0 + 12.0 * y as f32, 500.0 + 5.0 * (x + y) as f32]; + let (tx_n, ty_n) = (w.div_ceil(tw), h.div_ceil(tw)); + let mut file = vec![0u8; 16]; + let (mut offsets, mut counts) = (vec![0u64; tx_n * ty_n], vec![0u64; tx_n * ty_n]); + for ty in 0..ty_n { + for tx in 0..tx_n { + // a tile overhanging the image repeats its edge + let px = |x: usize, y: usize| rgb_at((tx * tw + x).min(w - 1), (ty * tw + y).min(h - 1)); + let luma: Vec = + (0..tw * tw).map(|i| to_ycc(px(i % tw, i / tw), gains)[0].round() as u16).collect(); + let (cw, ch) = (tw / hs, tw / vs); + let chroma = |k: usize| -> Vec { (0..cw * ch).map(|i| to_ycc(px(i % cw * hs, i / cw * vs), gains)[k].round() as u16).collect() }; + let enc = ljpeg::encode_planar(&[luma, chroma(1), chroma(2)], tw, tw, &[(hs, vs), (1, 1), (1, 1)], 16, 1, 0); + offsets[ty * tx_n + tx] = file.len() as u64; + counts[ty * tx_n + tx] = enc.len() as u64; + file.extend(enc); + } + } + let info = ImageInfo { + width: w as u32, + height: h as u32, + bits_per_sample: vec![15, 15, 15], + samples_per_pixel: 3, + compression: 7, + photometric: photometric::YCBCR, + planar: 1, + predictor: 1, + sample_format: 1, + new_subfile_type: 0, + layout: Layout::Tiles { tile_width: tw as u32, tile_height: tw as u32 }, + offsets, + byte_counts: counts, + }; + let (luma, chroma, sampling) = read_ycbcr_planes(&file, &info).unwrap(); + assert_eq!(sampling, (hs, vs)); + let out = ycbcr_to_camera_rgb(&luma, &chroma, w, h, sampling, [0.299, 0.587, 0.114], gains); + for y in 0..h { + for x in 0..w { + let want = rgb_at(x, y); + for c in 0..3 { + let got = out[(y * w + x) * 3 + c] as f32 - YCC_OUT_BLACK; + // rounding of the samples, amplified by the inverse transform; the last odd column / row + // has no chroma sample to its right / below and repeats the previous one + let edge = (hs == 2 && x == w - 1 && x % 2 == 1) || (vs == 2 && y == h - 1 && y % 2 == 1); + let tol = if edge { 40.0 } else { 3.0 }; + assert!((got - want[c]).abs() <= tol, "{hs}×{vs} ({x}, {y}) channel {c}: {got} vs {}", want[c]); + } + } + } + // a corrupt tile is skipped (left at zero luma) as long as one decodes; none decodable is an error + let mut bad = info.clone(); + bad.offsets[0] = file.len() as u64 + 10; + assert!(read_ycbcr_planes(&file, &bad).is_ok()); + bad.offsets.iter_mut().for_each(|o| *o = 3); + assert!(read_ycbcr_planes(&file, &bad).is_err()); + for n in [0, 20, file.len() / 2] { + let _ = read_ycbcr_planes(&file[..n], &info); + } + } + } + #[test] fn quad_tiles_place_cfa_cells() { // 14×12 mosaic in 8×8 tiles (the right and bottom tiles overhang), each tile a 4×4 frame of 2×2 cells,