WIP Issue #147: Sony M/S-size lossless ARW (subsampled YCbCr) — interrupted, unverified, do not land

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

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);

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<String> = 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<u16>)> = 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();
+}

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<u8>,
+    /// Horizontal and vertical sampling factors per component (T.81 `Hi`, `Vi`).
+    sampling: Vec<(usize, usize)>,
     tables: [Option<Huffman>; 4],
     table_for: Vec<usize>,
     predictor: u8,
@@ -228,6 +232,7 @@ fn parse_header(d: &[u8]) -> Result<Header, RawError> {
         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<Header, RawError> {
                 }
                 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::<usize>() > 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<Vec<(usize, usize)>, 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<Frame, RawError> {
     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<Frame, RawError> {
     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<u16>,
+}
+
+/// 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<Plane>,
+}
+
+/// 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<PlanarFrame, RawError> {
+    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::<Result<_, RawError>>()?;
+    let mut planes: Vec<Vec<u16>> = 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<u16>],
+    w: usize,
+    h: usize,
+    sampling: &[(usize, usize)],
+    precision: u8,
+    predictor: u8,
+    restart: usize,
+) -> Vec<u8> {
+    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<u16>)> = 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<u8>)> = (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<Vec<u16>> {
+        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::<Vec<_>>(), 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

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<Vec<u16>> {
     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<RawImage> {
+    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<u16>, [Vec<u16>; 2], (usize, usize));
+fn read_ycbcr_planes(bytes: &[u8], info: &ImageInfo) -> Result<YccPlanes> {
+    let (w, h) = (info.width as usize, info.height as usize);
+    let chunks = info.chunks(bytes.len() as u64);
+    let decoded: Vec<Result<(Chunk, ljpeg::PlanarFrame)>> = 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::<Vec<_>>()
+                )));
+            }
+            Ok((*c, f))
+        })
+        .collect();
+    let mut sampling = None;
+    let mut luma = vec![0u16; w * h];
+    let mut chroma: [Vec<u16>; 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<u16>; 2], w: usize, h: usize, (hs, vs): (usize, usize), coeffs: [f32; 3], gains: [f32; 3]) -> Vec<u16> {
+    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<u16>| {
+                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<RawImage> {
     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<RawImage> {
         _ => 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<u16> =
+                        (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<u16> { (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,