diff --git a/Cargo.lock b/Cargo.lock
index d6b9c020..e8f70c57 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -2,6 +2,12 @@
 # It is not intended for manual editing.
 version = 4
 
+[[package]]
+name = "adler2"
+version = "2.0.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "320119579fcad9c21884f5c4861d16174d0e06250625266f50fe6898340abefa"
+
 [[package]]
 name = "autocfg"
 version = "1.5.1"
@@ -166,6 +172,18 @@ dependencies = [
 [[package]]
 name = "lightcraft-raw"
 version = "0.1.0"
+dependencies = [
+ "lightcraft-color",
+ "lightcraft-geom",
+ "lightcraft-meta",
+ "lightcraft-raster",
+ "lightcraft-tiff",
+ "miniz_oxide",
+ "proptest",
+ "rayon",
+ "serde",
+ "thiserror",
+]
 
 [[package]]
 name = "lightcraft-tiff"
@@ -188,6 +206,15 @@ version = "2.8.3"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "cf8baf1c55e62ffcace7a9f06f4bd9cd3f0c4beb022d3b367256b91b87513d98"
 
+[[package]]
+name = "miniz_oxide"
+version = "0.8.9"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "1fa76a2c86f704bdb222d66965fb3d63269ce38518b83cb0575fca855ebb6316"
+dependencies = [
+ "adler2",
+]
+
 [[package]]
 name = "num-traits"
 version = "0.2.19"
diff --git a/crates/raw/Cargo.toml b/crates/raw/Cargo.toml
index 8cafc30f..512258bf 100644
--- a/crates/raw/Cargo.toml
+++ b/crates/raw/Cargo.toml
@@ -9,3 +9,15 @@ description = "LightCraft RAW decoding, demosaicing and DNG colour model (L1)"
 workspace = true
 
 [dependencies]
+lightcraft-tiff = { workspace = true }
+lightcraft-meta = { workspace = true }
+lightcraft-raster = { workspace = true }
+lightcraft-color = { workspace = true }
+lightcraft-geom = { workspace = true }
+serde = { workspace = true }
+thiserror = { workspace = true }
+rayon = { workspace = true }
+miniz_oxide = { workspace = true }
+
+[dev-dependencies]
+proptest = { workspace = true }
diff --git a/crates/raw/src/color.rs b/crates/raw/src/color.rs
new file mode 100644
index 00000000..5180745f
--- /dev/null
+++ b/crates/raw/src/color.rs
@@ -0,0 +1,377 @@
+//! DNG colour model (Adobe DNG Specification 1.7, chapter 6 "Mapping Camera Color Space to CIE XYZ Space").
+//!
+//! - `ColorMatrix1/2` (XYZ → reference camera) and `CameraCalibration1/2` / `ForwardMatrix1/2` are interpolated
+//!   between the two calibration illuminants by inverse correlated colour temperature of the white point.
+//! - `XYZtoCamera = AnalogBalance · CameraCalibration · ColorMatrix`.
+//! - White balance multipliers are the reciprocal of the camera neutral for the chosen white.
+//! - Camera → XYZ(D50): through `ForwardMatrix` when present, else the inverse of `XYZtoCamera` followed by
+//!   Bradford adaptation from the white to D50.
+//! - The as-shot white comes from `AsShotWhiteXY`, or `AsShotNeutral` converted to xy by the spec's iteration.
+//!
+//! Cameras without a colour matrix (non-DNG files until our own calibration DB exists) use a generic
+//! "camera RGB ≈ linear sRGB" model and are flagged with `matrix_is_fallback`.
+
+use crate::{ColorData, Mat3, RawImage};
+use lightcraft_color::{D50, D65, REC2020, SRGB, Xy, bradford, cct};
+use serde::{Deserialize, Serialize};
+
+/// Correlated colour temperature (K) of an Exif `LightSource` / DNG `CalibrationIlluminant` code.
+pub fn illuminant_temperature(code: u16) -> Option<f64> {
+    Some(match code {
+        1 | 4 | 9 => 5500.0,
+        2 | 14 => 4150.0,
+        3 => 2850.0,
+        10 => 6500.0,
+        11 => 7500.0,
+        12 => 6430.0,
+        13 => 5000.0,
+        15 => 3450.0,
+        16 => 2940.0,
+        17 => 2856.0,
+        18 => 4874.0,
+        19 => 6774.0,
+        20 => 5503.0,
+        21 => 6504.0,
+        22 => 7504.0,
+        23 => 5003.0,
+        24 => 3200.0,
+        _ => return None,
+    })
+}
+
+fn lerp(a: &Mat3, b: &Mat3, g: f64) -> Mat3 {
+    let mut r = [[0.0; 3]; 3];
+    for i in 0..3 {
+        for j in 0..3 {
+            r[i][j] = g * a.0[i][j] + (1.0 - g) * b.0[i][j];
+        }
+    }
+    Mat3(r)
+}
+
+/// Weight of calibration 1 for white `xy` (1 = use matrix 1 only).
+pub fn illuminant_weight(color: &ColorData, white: Xy) -> f64 {
+    let t = cct::xy_to_temp_tint(white).0;
+    let (t1, t2) = (illuminant_temperature(color.illuminant[0]).unwrap_or(2856.0), illuminant_temperature(color.illuminant[1]).unwrap_or(6504.0));
+    if (t1 - t2).abs() < 1.0 {
+        return 1.0;
+    }
+    let g = (1.0 / t - 1.0 / t2) / (1.0 / t1 - 1.0 / t2);
+    g.clamp(0.0, 1.0)
+}
+
+fn interp(pair: &[Option<Mat3>; 2], g: f64) -> Option<Mat3> {
+    match pair {
+        [Some(a), Some(b)] => Some(lerp(a, b, g)),
+        [Some(a), None] | [None, Some(a)] => Some(*a),
+        [None, None] => None,
+    }
+}
+
+/// `AnalogBalance · CameraCalibration` for white `xy`.
+fn ab_cc(color: &ColorData, g: f64) -> Mat3 {
+    let ab = color.analog_balance.map(|a| Mat3::diag(a[0], a[1], a[2])).unwrap_or(Mat3::IDENTITY);
+    let cc = interp(&color.camera_calibration, g).unwrap_or(Mat3::IDENTITY);
+    ab.mul(&cc)
+}
+
+/// Whether the file carries a usable colour matrix.
+pub fn has_matrix(color: &ColorData) -> bool {
+    color.color_matrix.iter().any(|m| m.is_some())
+}
+
+/// `XYZtoCamera` for white `xy` (fallback: XYZ → linear sRGB).
+pub fn xyz_to_camera(color: &ColorData, white: Xy) -> Mat3 {
+    let g = illuminant_weight(color, white);
+    match interp(&color.color_matrix, g) {
+        Some(cm) => ab_cc(color, g).mul(&cm),
+        None => SRGB.from_xyz(),
+    }
+}
+
+/// Camera neutral (raw camera values of a neutral surface under white `xy`), normalised to max 1.
+pub fn camera_neutral(color: &ColorData, white: Xy) -> [f64; 3] {
+    let n = xyz_to_camera(color, white).apply(white.to_xyz());
+    let mx = n.iter().cloned().fold(f64::MIN, f64::max);
+    if mx > 0.0 { n.map(|v| v / mx) } else { [1.0; 3] }
+}
+
+/// White-balance multipliers for white `xy` (reciprocal camera neutral, minimum 1).
+pub fn wb_multipliers(color: &ColorData, white: Xy) -> [f64; 3] {
+    let n = camera_neutral(color, white);
+    let m = n.map(|v| if v > 1e-9 { 1.0 / v } else { 1.0 });
+    let mn = m.iter().cloned().fold(f64::MAX, f64::min);
+    m.map(|v| v / mn)
+}
+
+/// Convert a camera neutral to a white xy (DNG spec iteration: xy → XYZtoCamera(xy) → xy until stable).
+pub fn neutral_to_xy(color: &ColorData, neutral: [f64; 3]) -> Xy {
+    let mut xy = D50;
+    for _ in 0..50 {
+        let Some(inv) = xyz_to_camera(color, xy).inverse() else { return xy };
+        let xyz = inv.apply(neutral);
+        let next = Xy::from_xyz(xyz);
+        if !(next.x.is_finite() && next.y.is_finite()) || next.y <= 0.0 {
+            return xy;
+        }
+        let next = Xy::new(next.x.clamp(0.1, 0.7), next.y.clamp(0.1, 0.7));
+        let done = (next.x - xy.x).abs() < 1e-9 && (next.y - xy.y).abs() < 1e-9;
+        xy = next;
+        if done {
+            break;
+        }
+    }
+    xy
+}
+
+/// The as-shot white point: `AsShotWhiteXY`, `AsShotNeutral`, vendor multipliers, else D65.
+pub fn as_shot_white_xy(raw: &RawImage) -> Xy {
+    let c = &raw.color;
+    if let Some(xy) = c.as_shot_white_xy {
+        return xy;
+    }
+    if let Some(n) = c.as_shot_neutral {
+        return neutral_to_xy(c, n);
+    }
+    if let Some(m) = raw.wb_multipliers.filter(|m| m.iter().all(|v| *v > 0.0 && v.is_finite())) {
+        return neutral_to_xy(c, [1.0 / m[0] as f64, 1.0 / m[1] as f64, 1.0 / m[2] as f64]);
+    }
+    D65
+}
+
+/// Camera (raw, not white balanced) → XYZ relative to D50, for scene white `xy`.
+pub fn camera_to_xyz_d50(color: &ColorData, white: Xy) -> Mat3 {
+    let g = illuminant_weight(color, white);
+    if has_matrix(color)
+        && let Some(fm) = interp(&color.forward_matrix, g)
+    {
+        let abcc = ab_cc(color, g);
+        if let Some(abcc_inv) = abcc.inverse() {
+            let n = camera_neutral(color, white);
+            let rn = abcc_inv.apply(n);
+            if rn.iter().all(|v| *v > 1e-9) {
+                return fm.mul(&Mat3::diag(1.0 / rn[0], 1.0 / rn[1], 1.0 / rn[2])).mul(&abcc_inv);
+            }
+        }
+    }
+    let cam_to_xyz = xyz_to_camera(color, white).inverse().unwrap_or(Mat3::IDENTITY);
+    bradford(white, D50).mul(&cam_to_xyz)
+}
+
+/// Everything the pipeline needs to go from raw camera RGB to the working space.
+#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]
+pub struct CameraTransform {
+    /// White-balanced camera RGB → linear Rec.2020 D65; maps (1, 1, 1) to (1, 1, 1).
+    pub matrix: Mat3,
+    /// Per-channel white-balance multipliers applied before `matrix` (minimum 1).
+    pub wb: [f32; 3],
+    pub white_xy: Xy,
+    /// No colour matrix was available: a generic camera ≈ sRGB model was used.
+    pub matrix_is_fallback: bool,
+    pub baseline_exposure: f64,
+}
+
+/// Camera transform for white `wb_xy`: `rec2020 = matrix · (wb ⊙ camera)`.
+pub fn camera_transform(raw: &RawImage, wb_xy: Xy) -> CameraTransform {
+    let color = &raw.color;
+    let fallback = !has_matrix(color);
+    let wb = wb_multipliers(color, wb_xy);
+    let to_d50 = camera_to_xyz_d50(color, wb_xy);
+    let m = REC2020.from_xyz().mul(&bradford(D50, D65)).mul(&to_d50).mul(&Mat3::diag(1.0 / wb[0], 1.0 / wb[1], 1.0 / wb[2]));
+    let white = m.apply([1.0; 3]);
+    let k = (white[0] + white[1] + white[2]) / 3.0;
+    let matrix = if k.abs() > 1e-12 { Mat3(m.0.map(|r| r.map(|v| v / k))) } else { m };
+    CameraTransform { matrix, wb: wb.map(|v| v as f32), white_xy: wb_xy, matrix_is_fallback: fallback, baseline_exposure: color.baseline_exposure }
+}
+
+/// `(matrix, multipliers)` with `rec2020 = matrix · (multipliers ⊙ camera)`, linear Rec.2020 D65.
+pub fn camera_to_rec2020(raw: &RawImage, wb_xy: Xy) -> (Mat3, [f32; 3]) {
+    let t = camera_transform(raw, wb_xy);
+    (t.matrix, t.wb)
+}
+
+/// Grey-world white balance estimate (multipliers, min 1) from camera RGB, ignoring clipped and very dark pixels.
+pub fn grey_world(img: &lightcraft_raster::Rgb32f) -> [f32; 3] {
+    let mut s = [0f64; 3];
+    let step = (img.data.len() / 200_000).max(1);
+    for p in img.data.iter().step_by(step) {
+        if p.iter().all(|&v| v > 0.02 && v < 0.95) {
+            for c in 0..3 {
+                s[c] += p[c] as f64;
+            }
+        }
+    }
+    if s.iter().any(|&v| v <= 0.0) {
+        return [1.0; 3];
+    }
+    let m = [s[1] / s[0], 1.0, s[1] / s[2]];
+    let mn = m.iter().cloned().fold(f64::MAX, f64::min);
+    m.map(|v| (v / mn) as f32)
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    /// A plausible camera: XYZ → camera matrix (rows sum to something sensible), invertible.
+    fn cam_a() -> Mat3 {
+        Mat3([[0.9, 0.2, -0.15], [-0.3, 1.25, 0.08], [0.02, -0.12, 0.85]])
+    }
+    fn cam_b() -> Mat3 {
+        Mat3([[0.7, 0.3, -0.1], [-0.35, 1.3, 0.1], [0.05, -0.2, 1.0]])
+    }
+
+    fn color_single(m: Mat3, ill: u16) -> ColorData {
+        ColorData { illuminant: [ill, 0], color_matrix: [Some(m), None], ..Default::default() }
+    }
+
+    fn raw_with(color: ColorData) -> RawImage {
+        RawImage {
+            format: crate::RawFormat::Dng,
+            width: 1,
+            height: 1,
+            cpp: 1,
+            data: crate::RawData::U16(vec![0]),
+            cfa: None,
+            bits: 16,
+            black: Default::default(),
+            white: vec![65535.0],
+            active_area: crate::Rect::new(0, 0, 1, 1),
+            crop: crate::Rect::new(0, 0, 1, 1),
+            orientation: Default::default(),
+            color,
+            wb_multipliers: None,
+            linearized: false,
+            opcodes: Default::default(),
+            metadata: Default::default(),
+        }
+    }
+
+    fn close(a: [f64; 3], b: [f64; 3], tol: f64) -> bool {
+        (0..3).all(|i| (a[i] - b[i]).abs() < tol)
+    }
+
+    #[test]
+    fn neutral_maps_to_white_and_primaries_are_recovered() {
+        let raw = raw_with(color_single(cam_a(), 21));
+        let t = camera_transform(&raw, D65);
+        assert!(!t.matrix_is_fallback);
+        // a surface with XYZ = Rec.2020 red primary under D65, seen by the camera
+        let red_xyz = REC2020.to_xyz().apply([1.0, 0.0, 0.0]);
+        let cam = cam_a().apply(red_xyz);
+        let wbd = [cam[0] * t.wb[0] as f64, cam[1] * t.wb[1] as f64, cam[2] * t.wb[2] as f64];
+        let out = t.matrix.apply(wbd);
+        let k = out[0];
+        assert!(close([out[0] / k, out[1] / k, out[2] / k], [1.0, 0.0, 0.0], 1e-6), "{out:?}");
+        // white
+        let n = cam_a().apply(D65.to_xyz());
+        let wbd = [n[0] * t.wb[0] as f64, n[1] * t.wb[1] as f64, n[2] * t.wb[2] as f64];
+        let out = t.matrix.apply(wbd);
+        assert!((out[0] - out[1]).abs() < 1e-6 && (out[1] - out[2]).abs() < 1e-6, "{out:?}");
+        assert!(close(t.matrix.apply([1.0; 3]), [1.0; 3], 1e-12));
+        assert!(t.wb.iter().cloned().fold(f32::MAX, f32::min) == 1.0);
+    }
+
+    #[test]
+    fn dual_illuminant_interpolation() {
+        let color = ColorData { illuminant: [17, 21], color_matrix: [Some(cam_a()), Some(cam_b())], ..Default::default() };
+        assert!((illuminant_weight(&color, cct::temp_tint_to_xy(2856.0, 0.0)) - 1.0).abs() < 0.01);
+        assert!(illuminant_weight(&color, cct::temp_tint_to_xy(6504.0, 0.0)) < 0.01);
+        assert_eq!(illuminant_weight(&color, cct::temp_tint_to_xy(10000.0, 0.0)), 0.0);
+        let mid = cct::temp_tint_to_xy(4000.0, 0.0);
+        let g = illuminant_weight(&color, mid);
+        let expect = (1.0 / 4000.0 - 1.0 / 6504.0) / (1.0 / 2856.0 - 1.0 / 6504.0);
+        assert!((g - expect).abs() < 0.02, "{g} vs {expect}");
+        let m = xyz_to_camera(&color, mid);
+        let l = lerp(&cam_a(), &cam_b(), g);
+        assert!(close(m.0[0], l.0[0], 1e-12));
+    }
+
+    #[test]
+    fn neutral_xy_roundtrip() {
+        let color = ColorData {
+            illuminant: [17, 21],
+            color_matrix: [Some(cam_a()), Some(cam_b())],
+            analog_balance: Some([1.0, 1.05, 0.98]),
+            ..Default::default()
+        };
+        for t in [2600.0, 3200.0, 4500.0, 5500.0, 7000.0, 9000.0] {
+            for tint in [-10.0, 0.0, 15.0] {
+                let xy = cct::temp_tint_to_xy(t, tint);
+                let n = camera_neutral(&color, xy);
+                let back = neutral_to_xy(&color, n);
+                assert!((back.x - xy.x).abs() < 1e-6 && (back.y - xy.y).abs() < 1e-6, "{t} {tint}: {back:?} vs {xy:?}");
+            }
+        }
+        let mut raw = raw_with(color.clone());
+        raw.color.as_shot_neutral = Some(camera_neutral(&color, cct::temp_tint_to_xy(3300.0, 5.0)));
+        let xy = as_shot_white_xy(&raw);
+        let (t, tint) = cct::xy_to_temp_tint(xy);
+        assert!((t - 3300.0).abs() < 5.0 && (tint - 5.0).abs() < 0.5, "{t} {tint}");
+        raw.color.as_shot_white_xy = Some(D50);
+        assert_eq!(as_shot_white_xy(&raw), D50);
+    }
+
+    #[test]
+    fn forward_matrix_path_agrees_with_color_matrix_path() {
+        let a = cam_a();
+        let n = a.apply(D65.to_xyz());
+        let fm = bradford(D65, D50).mul(&a.inverse().unwrap()).mul(&Mat3::diag(n[0], n[1], n[2]));
+        // DNG forward matrices map (1,1,1) to D50 white (with Y=1): normalise like a real file would
+        let w = fm.apply([1.0; 3]);
+        let s = D50.to_xyz()[1] / w[1];
+        let fm = Mat3(fm.0.map(|r| r.map(|v| v * s)));
+        let cm_only = camera_transform(&raw_with(color_single(a, 21)), D65);
+        let mut with_fm = color_single(a, 21);
+        with_fm.forward_matrix = [Some(fm), None];
+        let fm_t = camera_transform(&raw_with(with_fm), D65);
+        for i in 0..3 {
+            assert!(close(cm_only.matrix.0[i], fm_t.matrix.0[i], 1e-9), "{:?} vs {:?}", cm_only.matrix, fm_t.matrix);
+        }
+        assert_eq!(cm_only.wb, fm_t.wb);
+    }
+
+    #[test]
+    fn camera_calibration_and_analog_balance_are_respected() {
+        // AB · CC scales channels: the neutral must scale accordingly
+        let mut c = color_single(cam_a(), 21);
+        let base = camera_neutral(&c, D65);
+        c.camera_calibration = [Some(Mat3::diag(1.0, 0.5, 1.0)), None];
+        let n = camera_neutral(&c, D65);
+        assert!((n[1] / n[0] - 0.5 * base[1] / base[0]).abs() < 1e-9);
+    }
+
+    #[test]
+    fn fallback_without_matrix() {
+        let mut raw = raw_with(ColorData::default());
+        raw.wb_multipliers = Some([2.0, 1.0, 1.5]);
+        let xy = as_shot_white_xy(&raw);
+        let t = camera_transform(&raw, xy);
+        assert!(t.matrix_is_fallback);
+        assert!((t.wb[0] / t.wb[1] - 2.0).abs() < 1e-3 && (t.wb[2] / t.wb[1] - 1.5).abs() < 1e-3, "{:?}", t.wb);
+        let t = camera_transform(&raw, D65);
+        let expect = SRGB.to_space(&REC2020);
+        for i in 0..3 {
+            assert!(close(t.matrix.0[i], expect.0[i], 1e-6));
+        }
+        assert_eq!(camera_to_rec2020(&raw, D65).1, [1.0, 1.0, 1.0]);
+    }
+
+    #[test]
+    fn illuminant_codes() {
+        assert_eq!(illuminant_temperature(21), Some(6504.0));
+        assert_eq!(illuminant_temperature(17), Some(2856.0));
+        assert_eq!(illuminant_temperature(0), None);
+        assert_eq!(illuminant_temperature(255), None);
+    }
+
+    #[test]
+    fn grey_world_balances() {
+        let img = lightcraft_raster::Rgb32f::from_fn(20, 20, |x, _| {
+            let v = 0.2 + 0.02 * x as f32;
+            [v * 0.5, v, v * 0.8]
+        });
+        let m = grey_world(&img);
+        assert!((m[0] - 2.0).abs() < 1e-3 && (m[1] - 1.0).abs() < 1e-6 && (m[2] - 1.25).abs() < 1e-3, "{m:?}");
+    }
+}
diff --git a/crates/raw/src/demosaic/ahd.rs b/crates/raw/src/demosaic/ahd.rs
new file mode 100644
index 00000000..aaf50bc6
--- /dev/null
+++ b/crates/raw/src/demosaic/ahd.rs
@@ -0,0 +1,165 @@
+//! Adaptive Homogeneity-Directed demosaicing (K. Hirakawa, T. W. Parks, "Adaptive homogeneity-directed
+//! demosaicing algorithm", IEEE Trans. Image Processing 14(3), 2005), implemented from the paper:
+//!
+//! 1. Interpolate green horizontally and vertically (5-tap filter [−¼ ½ ½ ½ −¼] across the missing site,
+//!    clamped to the two adjacent greens), giving two candidate green planes.
+//! 2. Complete each candidate with red/blue by colour-difference (C − G) interpolation.
+//! 3. Convert both candidates to CIELab and measure, per pixel, how many of its 4 neighbours lie within the
+//!    adaptive luminance / chrominance tolerances (ε_L, ε_C from the paper) — the homogeneity.
+//! 4. Sum homogeneity over a 3×3 window and pick the more homogeneous candidate (average on ties).
+//!
+//! Work is done in overlapping tiles (bounded memory, rayon-parallel).
+
+use super::Mosaic;
+use crate::Rgb32f;
+use rayon::prelude::*;
+
+const TILE: usize = 128;
+const MARGIN: usize = 6;
+
+pub(crate) fn ahd(m: &Mosaic) -> Rgb32f {
+    let (w, h) = (m.w, m.h);
+    let tiles: Vec<(usize, usize)> = (0..h.div_ceil(TILE)).flat_map(|ty| (0..w.div_ceil(TILE)).map(move |tx| (tx * TILE, ty * TILE))).collect();
+    let results: Vec<(usize, usize, usize, usize, Vec<[f32; 3]>)> = tiles
+        .par_iter()
+        .map(|&(x0, y0)| {
+            let (tw, th) = (TILE.min(w - x0), TILE.min(h - y0));
+            (x0, y0, tw, th, tile(m, x0, y0, tw, th))
+        })
+        .collect();
+    let mut out = Rgb32f::new(w, h);
+    for (x0, y0, tw, th, px) in results {
+        for y in 0..th {
+            out.data[(y0 + y) * w + x0..(y0 + y) * w + x0 + tw].copy_from_slice(&px[y * tw..(y + 1) * tw]);
+        }
+    }
+    out
+}
+
+/// Linear camera RGB (treated as sRGB/D65 primaries for the homogeneity metric only) → CIELab.
+#[inline]
+fn lab(p: [f32; 3]) -> [f32; 3] {
+    let x = (0.4124 * p[0] + 0.3576 * p[1] + 0.1805 * p[2]) / 0.95047;
+    let y = 0.2126 * p[0] + 0.7152 * p[1] + 0.0722 * p[2];
+    let z = (0.0193 * p[0] + 0.1192 * p[1] + 0.9505 * p[2]) / 1.08883;
+    #[inline]
+    fn f(t: f32) -> f32 {
+        if t > 0.008856 { t.cbrt() } else { 7.787 * t + 16.0 / 116.0 }
+    }
+    let (fx, fy, fz) = (f(x), f(y), f(z));
+    [116.0 * fy - 16.0, 500.0 * (fx - fy), 200.0 * (fy - fz)]
+}
+
+fn tile(m: &Mosaic, x0: usize, y0: usize, tw: usize, th: usize) -> Vec<[f32; 3]> {
+    let (bw, bh) = (tw + 2 * MARGIN, th + 2 * MARGIN);
+    let n = bw * bh;
+    let (ox, oy) = (x0 as isize - MARGIN as isize, y0 as isize - MARGIN as isize);
+    let mut mos = vec![0f32; n];
+    let mut col = vec![0u8; n];
+    for y in 0..bh {
+        for x in 0..bw {
+            mos[y * bw + x] = m.get(ox + x as isize, oy + y as isize);
+            col[y * bw + x] = m.color(ox + x as isize, oy + y as isize);
+        }
+    }
+    let cl = |x: isize, y: isize| -> usize { y.clamp(0, bh as isize - 1) as usize * bw + x.clamp(0, bw as isize - 1) as usize };
+    // 1. directional greens
+    let mut gh = vec![0f32; n];
+    let mut gv = vec![0f32; n];
+    for y in 0..bh as isize {
+        for x in 0..bw as isize {
+            let i = cl(x, y);
+            let c = mos[i];
+            if col[i] == 1 {
+                gh[i] = c;
+                gv[i] = c;
+                continue;
+            }
+            let (l, r) = (mos[cl(x - 1, y)], mos[cl(x + 1, y)]);
+            let est = (l + r) * 0.5 + (2.0 * c - mos[cl(x - 2, y)] - mos[cl(x + 2, y)]) * 0.25;
+            gh[i] = est.clamp(l.min(r), l.max(r));
+            let (u, d) = (mos[cl(x, y - 1)], mos[cl(x, y + 1)]);
+            let est = (u + d) * 0.5 + (2.0 * c - mos[cl(x, y - 2)] - mos[cl(x, y + 2)]) * 0.25;
+            gv[i] = est.clamp(u.min(d), u.max(d));
+        }
+    }
+    // 2. full candidates
+    let complete = |g: &[f32]| -> Vec<[f32; 3]> {
+        let mut rgb = vec![[0f32; 3]; n];
+        for y in 0..bh as isize {
+            for x in 0..bw as isize {
+                let i = cl(x, y);
+                let own = col[i] as usize;
+                let gvv = g[i];
+                let d = |dx: isize, dy: isize| {
+                    let j = cl(x + dx, y + dy);
+                    mos[j] - g[j]
+                };
+                let p = &mut rgb[i];
+                p[1] = gvv;
+                if own == 1 {
+                    let ch = col[cl(x + 1, y)] as usize;
+                    let cv = col[cl(x, y + 1)] as usize;
+                    p[ch] = gvv + (d(-1, 0) + d(1, 0)) * 0.5;
+                    p[cv] = gvv + (d(0, -1) + d(0, 1)) * 0.5;
+                } else {
+                    p[own] = mos[i];
+                    p[2 - own] = gvv + (d(-1, -1) + d(1, -1) + d(-1, 1) + d(1, 1)) * 0.25;
+                }
+            }
+        }
+        rgb
+    };
+    let rh = complete(&gh);
+    let rv = complete(&gv);
+    // 3. Lab + homogeneity
+    let lh: Vec<[f32; 3]> = rh.iter().map(|&p| lab(p)).collect();
+    let lv: Vec<[f32; 3]> = rv.iter().map(|&p| lab(p)).collect();
+    let dl = |a: &[f32; 3], b: &[f32; 3]| (a[0] - b[0]).abs();
+    let dc = |a: &[f32; 3], b: &[f32; 3]| (a[1] - b[1]).powi(2) + (a[2] - b[2]).powi(2);
+    let mut hh = vec![0u8; n];
+    let mut hv = vec![0u8; n];
+    for y in 1..bh as isize - 1 {
+        for x in 1..bw as isize - 1 {
+            let i = cl(x, y);
+            let (l, r, u, d) = (cl(x - 1, y), cl(x + 1, y), cl(x, y - 1), cl(x, y + 1));
+            let eps_l = dl(&lh[i], &lh[l]).max(dl(&lh[i], &lh[r])).min(dl(&lv[i], &lv[u]).max(dl(&lv[i], &lv[d])));
+            let eps_c = dc(&lh[i], &lh[l]).max(dc(&lh[i], &lh[r])).min(dc(&lv[i], &lv[u]).max(dc(&lv[i], &lv[d])));
+            let mut ch = 0;
+            let mut cv = 0;
+            for j in [l, r, u, d] {
+                if dl(&lh[i], &lh[j]) <= eps_l && dc(&lh[i], &lh[j]) <= eps_c {
+                    ch += 1;
+                }
+                if dl(&lv[i], &lv[j]) <= eps_l && dc(&lv[i], &lv[j]) <= eps_c {
+                    cv += 1;
+                }
+            }
+            hh[i] = ch;
+            hv[i] = cv;
+        }
+    }
+    // 4. choose
+    let mut out = Vec::with_capacity(tw * th);
+    for y in MARGIN..MARGIN + th {
+        for x in MARGIN..MARGIN + tw {
+            let (mut sh, mut sv) = (0u32, 0u32);
+            for dy in -1isize..=1 {
+                for dx in -1isize..=1 {
+                    let j = cl(x as isize + dx, y as isize + dy);
+                    sh += hh[j] as u32;
+                    sv += hv[j] as u32;
+                }
+            }
+            let i = y * bw + x;
+            out.push(if sh > sv {
+                rh[i]
+            } else if sv > sh {
+                rv[i]
+            } else {
+                [(rh[i][0] + rv[i][0]) * 0.5, (rh[i][1] + rv[i][1]) * 0.5, (rh[i][2] + rv[i][2]) * 0.5]
+            });
+        }
+    }
+    out
+}
diff --git a/crates/raw/src/demosaic/bilinear.rs b/crates/raw/src/demosaic/bilinear.rs
new file mode 100644
index 00000000..cc614507
--- /dev/null
+++ b/crates/raw/src/demosaic/bilinear.rs
@@ -0,0 +1,78 @@
+//! Bilinear demosaicing for any CFA pattern: each missing colour is the mean of the same-colour samples in the
+//! 3×3 neighbourhood (5×5 when the 3×3 has none, as happens in some non-Bayer layouts).
+
+use super::Mosaic;
+use crate::Rgb32f;
+use lightcraft_raster::par_rows;
+
+/// Offsets (dx, dy) of same-colour neighbours for each pattern position and colour.
+pub(crate) fn neighbour_table(m: &Mosaic) -> Vec<[Vec<(isize, isize)>; 3]> {
+    let cfa = m.cfa;
+    let mut table = Vec::with_capacity(cfa.width * cfa.height);
+    for py in 0..cfa.height {
+        for px in 0..cfa.width {
+            let mut per: [Vec<(isize, isize)>; 3] = Default::default();
+            for (c, list) in per.iter_mut().enumerate() {
+                for r in [1isize, 2] {
+                    for dy in -r..=r {
+                        for dx in -r..=r {
+                            let (x, y) = (px as isize + dx + 6 * cfa.width as isize, py as isize + dy + 6 * cfa.height as isize);
+                            if cfa.color_at(x as usize, y as usize) as usize == c {
+                                list.push((dx, dy));
+                            }
+                        }
+                    }
+                    if !list.is_empty() {
+                        break;
+                    }
+                }
+            }
+            table.push(per);
+        }
+    }
+    table
+}
+
+pub(crate) fn bilinear(m: &Mosaic) -> Rgb32f {
+    let (w, h) = (m.w, m.h);
+    let table = neighbour_table(m);
+    let cfa = m.cfa;
+    let mut out = Rgb32f::new(w, h);
+    par_rows(&mut out.data, w, |y, row| {
+        for (x, px) in row.iter_mut().enumerate() {
+            let own = cfa.color_at(x, y) as usize;
+            let interior = x >= 2 && y >= 2 && x + 2 < w && y + 2 < h;
+            let entry = &table[(y % cfa.height) * cfa.width + (x % cfa.width)];
+            for c in 0..3 {
+                if c == own {
+                    px[c] = m.data[y * w + x];
+                    continue;
+                }
+                let (mut s, mut n) = (0.0f32, 0u32);
+                if interior {
+                    for &(dx, dy) in &entry[c] {
+                        s += m.data[(y as isize + dy) as usize * w + (x as isize + dx) as usize];
+                        n += 1;
+                    }
+                } else {
+                    for r in [1isize, 2] {
+                        for dy in -r..=r {
+                            for dx in -r..=r {
+                                let (nx, ny) = (x as isize + dx, y as isize + dy);
+                                if m.color(nx, ny) as usize == c {
+                                    s += m.get(nx, ny);
+                                    n += 1;
+                                }
+                            }
+                        }
+                        if n > 0 {
+                            break;
+                        }
+                    }
+                }
+                px[c] = if n > 0 { s / n as f32 } else { m.data[y * w + x] };
+            }
+        }
+    });
+    out
+}
diff --git a/crates/raw/src/demosaic/mod.rs b/crates/raw/src/demosaic/mod.rs
new file mode 100644
index 00000000..f00913c5
--- /dev/null
+++ b/crates/raw/src/demosaic/mod.rs
@@ -0,0 +1,253 @@
+//! Demosaicing: CFA mosaic → camera RGB ([`Rgb32f`]), values as normalised (white = 1.0).
+//!
+//! Methods (all from published descriptions, see the README):
+//! - [`Method::Bilinear`] — average of same-colour neighbours (drafts / thumbnails); works for any pattern.
+//! - [`Method::Ppg`] — Patterned Pixel Grouping (Chuan-kai Lin): gradient-selected green, colour-difference R/B.
+//! - [`Method::Ahd`] — Adaptive Homogeneity-Directed (Hirakawa & Parks, IEEE TIP 2005): horizontal/vertical
+//!   candidates chosen per pixel by CIELab homogeneity. The quality default for Bayer sensors.
+//! - Non-Bayer patterns (X-Trans, …) always use `xtrans::directional`, our own edge-weighted colour-difference
+//!   interpolation.
+//!
+//! Every method reproduces a constant-colour mosaic exactly.
+
+mod ahd;
+mod bilinear;
+mod ppg;
+mod xtrans;
+
+use crate::{Cfa, Normalized, Rgb32f};
+use serde::{Deserialize, Serialize};
+
+#[derive(Clone, Copy, Debug, Default, PartialEq, Eq, Hash, Serialize, Deserialize)]
+pub enum Method {
+    Bilinear,
+    Ppg,
+    #[default]
+    Ahd,
+}
+
+/// Mirror-reflected access to a single-plane mosaic (reflection keeps Bayer parity: −k → k, w−1+k → w−1−k).
+#[derive(Clone, Copy)]
+pub(crate) struct Mosaic<'a> {
+    pub w: usize,
+    pub h: usize,
+    pub data: &'a [f32],
+    pub cfa: &'a Cfa,
+}
+
+#[inline]
+pub(crate) fn reflect(i: isize, n: usize) -> usize {
+    if i >= 0 && (i as usize) < n {
+        return i as usize;
+    }
+    let n = n as isize;
+    if n == 1 {
+        return 0;
+    }
+    let period = 2 * (n - 1);
+    let mut i = i.rem_euclid(period);
+    if i >= n {
+        i = period - i;
+    }
+    i as usize
+}
+
+impl Mosaic<'_> {
+    #[inline]
+    pub fn get(&self, x: isize, y: isize) -> f32 {
+        self.data[reflect(y, self.h) * self.w + reflect(x, self.w)]
+    }
+    #[inline]
+    pub fn color(&self, x: isize, y: isize) -> u8 {
+        self.cfa.color_at(reflect(x, self.w), reflect(y, self.h))
+    }
+}
+
+/// Demosaic normalised raw data. `cpp == 3` data is returned as is; monochrome data is replicated to grey.
+pub fn demosaic(n: &Normalized, method: Method) -> Rgb32f {
+    let (w, h) = (n.width, n.height);
+    if w == 0 || h == 0 {
+        return Rgb32f::new(w, h);
+    }
+    match (&n.cfa, n.cpp) {
+        (_, 3) => Rgb32f { width: w, height: h, data: n.data.chunks_exact(3).map(|c| [c[0], c[1], c[2]]).collect() },
+        (_, cpp) if cpp != 1 => Rgb32f { width: w, height: h, data: n.data.chunks_exact(cpp).map(|c| [c[0], c[1], c[2]]).collect() },
+        (None, _) => Rgb32f { width: w, height: h, data: n.data.iter().map(|&v| [v; 3]).collect() },
+        (Some(cfa), _) => {
+            let m = Mosaic { w, h, data: &n.data, cfa };
+            if !cfa.is_bayer() {
+                return match method {
+                    Method::Bilinear => bilinear::bilinear(&m),
+                    _ => xtrans::directional(&m),
+                };
+            }
+            if w < 4 || h < 4 {
+                return bilinear::bilinear(&m);
+            }
+            match method {
+                Method::Bilinear => bilinear::bilinear(&m),
+                Method::Ppg => ppg::ppg(&m),
+                Method::Ahd => ahd::ahd(&m),
+            }
+        }
+    }
+}
+
+/// Build a mosaic from an RGB image by sampling each pixel's CFA colour (for tests and synthetic data).
+pub fn mosaic_from_rgb(img: &Rgb32f, cfa: &Cfa) -> Normalized {
+    let data = (0..img.height).flat_map(|y| (0..img.width).map(move |x| (x, y))).map(|(x, y)| img.get(x, y)[cfa.color_at(x, y) as usize]).collect();
+    Normalized { width: img.width, height: img.height, cpp: 1, data, cfa: Some(cfa.clone()) }
+}
+
+/// Peak signal-to-noise ratio (dB, peak 1.0) between two images over an interior region (`border` pixels excluded).
+pub fn psnr(a: &Rgb32f, b: &Rgb32f, border: usize) -> f64 {
+    let mut se = 0.0f64;
+    let mut n = 0usize;
+    for y in border..a.height.saturating_sub(border) {
+        for x in border..a.width.saturating_sub(border) {
+            let (p, q) = (a.get(x, y), b.get(x, y));
+            for c in 0..3 {
+                se += ((p[c] - q[c]) as f64).powi(2);
+                n += 1;
+            }
+        }
+    }
+    if n == 0 || se == 0.0 {
+        return f64::INFINITY;
+    }
+    10.0 * (1.0 / (se / n as f64)).log10()
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    fn smooth_scene(w: usize, h: usize) -> Rgb32f {
+        Rgb32f::from_fn(w, h, |x, y| {
+            let (fx, fy) = (x as f32 / w as f32, y as f32 / h as f32);
+            [0.2 + 0.6 * fx * fy, 0.3 + 0.4 * (fx * 3.0).sin().abs() * 0.5 + 0.2 * fy, 0.1 + 0.5 * (1.0 - fx) * (0.5 + 0.5 * (fy * 5.0).cos())]
+        })
+    }
+
+    /// A photo-like scene: sharp luminance structure (disc, stripes, diagonal edge, fine texture) modulating a
+    /// smoothly varying colour — edges are mostly luminance edges, as in natural images.
+    fn edgy_scene(w: usize, h: usize) -> Rgb32f {
+        Rgb32f::from_fn(w, h, |x, y| {
+            let (fx, fy) = (x as f32, y as f32);
+            let d = ((fx - w as f32 * 0.35).powi(2) + (fy - h as f32 * 0.4).powi(2)).sqrt();
+            let lum = if d < w as f32 * 0.18 {
+                0.85
+            } else if x > w * 2 / 3 && (y / 3) % 2 == 0 {
+                0.15
+            } else if fx + fy * 0.5 > w as f32 * 1.1 {
+                0.7
+            } else {
+                0.35 + 0.1 * ((fx * 0.7).sin() * (fy * 0.9).cos())
+            };
+            let tint = [0.9 + 0.1 * (fx / w as f32), 0.8, 0.55 + 0.3 * (fy / h as f32)];
+            [lum * tint[0], lum * tint[1], lum * tint[2]]
+        })
+    }
+
+    fn run(img: &Rgb32f, cfa: &Cfa, m: Method) -> Rgb32f {
+        demosaic(&mosaic_from_rgb(img, cfa), m)
+    }
+
+    #[test]
+    fn constant_colour_is_exact_for_all_methods_and_phases() {
+        let img = Rgb32f::filled(17, 13, [0.3, 0.55, 0.8]);
+        for pat in ["RGGB", "BGGR", "GRBG", "GBRG"] {
+            let cfa = Cfa::bayer(pat).unwrap();
+            for m in [Method::Bilinear, Method::Ppg, Method::Ahd] {
+                let out = run(&img, &cfa, m);
+                for p in &out.data {
+                    for c in 0..3 {
+                        assert!((p[c] - img.data[0][c]).abs() < 1e-5, "{pat} {m:?} {p:?}");
+                    }
+                }
+            }
+        }
+        let out = run(&Rgb32f::filled(19, 14, [0.3, 0.55, 0.8]), &Cfa::xtrans(), Method::Ahd);
+        assert!(out.data.iter().all(|p| (p[0] - 0.3).abs() < 1e-5 && (p[1] - 0.55).abs() < 1e-5 && (p[2] - 0.8).abs() < 1e-5));
+    }
+
+    #[test]
+    fn known_samples_are_preserved() {
+        let img = edgy_scene(40, 32);
+        let cfa = Cfa::bayer("GRBG").unwrap();
+        for m in [Method::Bilinear, Method::Ppg, Method::Ahd] {
+            let out = run(&img, &cfa, m);
+            for y in 0..32 {
+                for x in 0..40 {
+                    let c = cfa.color_at(x, y) as usize;
+                    assert!((out.get(x, y)[c] - img.get(x, y)[c]).abs() < 1e-6, "{m:?}");
+                }
+            }
+        }
+    }
+
+    #[test]
+    fn psnr_thresholds_and_ranking() {
+        let smooth = smooth_scene(96, 80);
+        let edgy = edgy_scene(96, 80);
+        let cfa = Cfa::bayer("RGGB").unwrap();
+        let p = |img: &Rgb32f, m| psnr(img, &run(img, &cfa, m), 4);
+        let (sb, sp, sa) = (p(&smooth, Method::Bilinear), p(&smooth, Method::Ppg), p(&smooth, Method::Ahd));
+        assert!(sb > 38.0 && sp > 38.0 && sa > 38.0, "smooth: bilinear {sb:.1} ppg {sp:.1} ahd {sa:.1}");
+        let (eb, ep, ea) = (p(&edgy, Method::Bilinear), p(&edgy, Method::Ppg), p(&edgy, Method::Ahd));
+        eprintln!("smooth: bilinear {sb:.1} ppg {sp:.1} ahd {sa:.1}; edgy: bilinear {eb:.1} ppg {ep:.1} ahd {ea:.1}");
+        assert!(eb > 20.0, "edgy bilinear {eb:.1}");
+        assert!(ep > eb + 6.0, "ppg {ep:.1} vs bilinear {eb:.1}");
+        assert!(ea > eb + 6.0, "ahd {ea:.1} vs bilinear {eb:.1}");
+        let xt = Cfa::xtrans();
+        let xp = psnr(&edgy, &run(&edgy, &xt, Method::Ahd), 6);
+        let xb = psnr(&edgy, &run(&edgy, &xt, Method::Bilinear), 6);
+        assert!(xp > 20.0 && xp >= xb - 0.5, "xtrans directional {xp:.1} bilinear {xb:.1}");
+    }
+
+    #[test]
+    fn tiny_and_degenerate_inputs() {
+        let cfa = Cfa::bayer("RGGB").unwrap();
+        for (w, h) in [(1, 1), (2, 2), (3, 5), (5, 3), (1, 7)] {
+            let img = Rgb32f::filled(w, h, [0.5; 3]);
+            for m in [Method::Bilinear, Method::Ppg, Method::Ahd] {
+                let out = run(&img, &cfa, m);
+                assert_eq!((out.width, out.height), (w, h));
+                assert!(out.data.iter().all(|p| p.iter().all(|v| v.is_finite())));
+            }
+        }
+        let n = Normalized { width: 2, height: 1, cpp: 1, data: vec![0.2, 0.4], cfa: None };
+        assert_eq!(demosaic(&n, Method::Ahd).data, vec![[0.2; 3], [0.4; 3]]);
+        let n = Normalized { width: 1, height: 1, cpp: 3, data: vec![0.1, 0.2, 0.3], cfa: None };
+        assert_eq!(demosaic(&n, Method::Ahd).data, vec![[0.1, 0.2, 0.3]]);
+    }
+
+    #[test]
+    fn reflect_indices() {
+        assert_eq!(reflect(-1, 5), 1);
+        assert_eq!(reflect(-2, 5), 2);
+        assert_eq!(reflect(5, 5), 3);
+        assert_eq!(reflect(6, 5), 2);
+        assert_eq!(reflect(3, 1), 0);
+        assert_eq!(reflect(-7, 2), 1);
+    }
+
+    #[test]
+    #[ignore]
+    fn throughput() {
+        let (w, h) = (6000, 4000);
+        let img = smooth_scene(w, h);
+        let cfa = Cfa::bayer("RGGB").unwrap();
+        let n = mosaic_from_rgb(&img, &cfa);
+        for m in [Method::Bilinear, Method::Ppg, Method::Ahd] {
+            let t = std::time::Instant::now();
+            let _ = demosaic(&n, m);
+            let s = t.elapsed().as_secs_f64();
+            println!("{m:?}: {:.1} MP/s ({:.0} ms for 24 MP)", (w * h) as f64 / 1e6 / s, s * 1e3);
+        }
+        let n = mosaic_from_rgb(&img, &Cfa::xtrans());
+        let t = std::time::Instant::now();
+        let _ = demosaic(&n, Method::Ahd);
+        println!("X-Trans directional: {:.1} MP/s", (w * h) as f64 / 1e6 / t.elapsed().as_secs_f64());
+    }
+}
diff --git a/crates/raw/src/demosaic/ppg.rs b/crates/raw/src/demosaic/ppg.rs
new file mode 100644
index 00000000..83779615
--- /dev/null
+++ b/crates/raw/src/demosaic/ppg.rs
@@ -0,0 +1,89 @@
+//! Patterned Pixel Grouping (after Chuan-kai Lin's published description of the method).
+//!
+//! 1. Green at red/blue sites: four directional gradients (same-colour difference ×2 plus the green difference
+//!    across the site); the estimate along the flattest direction is the nearer green weighted 3:1 with the far
+//!    one plus a same-colour Laplacian correction, clamped to the two greens of that direction.
+//! 2. Red/blue at green sites: colour-difference (C − G) average of the two same-colour neighbours.
+//! 3. Blue at red sites (and red at blue): colour-difference average along the diagonal with the smaller
+//!    gradient (both diagonals when equal).
+
+use super::Mosaic;
+use crate::Rgb32f;
+use lightcraft_raster::par_rows;
+
+/// Green plane with PPG estimates at non-green sites.
+pub(crate) fn green_plane(m: &Mosaic) -> Vec<f32> {
+    let (w, h) = (m.w, m.h);
+    let mut g = vec![0f32; w * h];
+    par_rows(&mut g, w, |y, row| {
+        let y = y as isize;
+        for (x, o) in row.iter_mut().enumerate() {
+            let x = x as isize;
+            let c = m.get(x, y);
+            if m.color(x, y) == 1 {
+                *o = c;
+                continue;
+            }
+            let (gn, gs, gw, ge) = (m.get(x, y - 1), m.get(x, y + 1), m.get(x - 1, y), m.get(x + 1, y));
+            let (cn, cs, cw, ce) = (m.get(x, y - 2), m.get(x, y + 2), m.get(x - 2, y), m.get(x + 2, y));
+            let dv = (gn - gs).abs();
+            let dh = (gw - ge).abs();
+            let grads = [(cn - c).abs() * 2.0 + dv, (cs - c).abs() * 2.0 + dv, (cw - c).abs() * 2.0 + dh, (ce - c).abs() * 2.0 + dh];
+            let mut best = 0;
+            for i in 1..4 {
+                if grads[i] < grads[best] {
+                    best = i;
+                }
+            }
+            let (near, far, cfar) = match best {
+                0 => (gn, gs, cn),
+                1 => (gs, gn, cs),
+                2 => (gw, ge, cw),
+                _ => (ge, gw, ce),
+            };
+            let est = (near * 3.0 + far + c - cfar) * 0.25;
+            *o = est.clamp(near.min(far), near.max(far));
+        }
+    });
+    g
+}
+
+pub(crate) fn ppg(m: &Mosaic) -> Rgb32f {
+    let (w, h) = (m.w, m.h);
+    let g = green_plane(m);
+    let gg = |x: isize, y: isize| g[super::reflect(y, h) * w + super::reflect(x, w)];
+    let mut out = Rgb32f::new(w, h);
+    par_rows(&mut out.data, w, |y, row| {
+        let y = y as isize;
+        for (x, px) in row.iter_mut().enumerate() {
+            let x = x as isize;
+            let own = m.color(x, y) as usize;
+            let gv = gg(x, y);
+            px[1] = gv;
+            if own == 1 {
+                // horizontal neighbours carry one colour, vertical the other
+                let ch = m.color(x + 1, y) as usize;
+                let cv = m.color(x, y + 1) as usize;
+                let hval = gv + ((m.get(x - 1, y) - gg(x - 1, y)) + (m.get(x + 1, y) - gg(x + 1, y))) * 0.5;
+                let vval = gv + ((m.get(x, y - 1) - gg(x, y - 1)) + (m.get(x, y + 1) - gg(x, y + 1))) * 0.5;
+                px[ch] = hval;
+                px[cv] = vval;
+            } else {
+                px[own] = m.get(x, y);
+                let other = 2 - own;
+                let d = |dx: isize, dy: isize| m.get(x + dx, y + dy) - gg(x + dx, y + dy);
+                let ne = (m.get(x + 1, y - 1) - m.get(x - 1, y + 1)).abs() + (gg(x + 1, y - 1) - gv).abs() + (gg(x - 1, y + 1) - gv).abs();
+                let nw = (m.get(x - 1, y - 1) - m.get(x + 1, y + 1)).abs() + (gg(x - 1, y - 1) - gv).abs() + (gg(x + 1, y + 1) - gv).abs();
+                let diff = if ne < nw {
+                    (d(1, -1) + d(-1, 1)) * 0.5
+                } else if nw < ne {
+                    (d(-1, -1) + d(1, 1)) * 0.5
+                } else {
+                    (d(1, -1) + d(-1, 1) + d(-1, -1) + d(1, 1)) * 0.25
+                };
+                px[other] = gv + diff;
+            }
+        }
+    });
+    out
+}
diff --git a/crates/raw/src/demosaic/xtrans.rs b/crates/raw/src/demosaic/xtrans.rs
new file mode 100644
index 00000000..c26f768b
--- /dev/null
+++ b/crates/raw/src/demosaic/xtrans.rs
@@ -0,0 +1,95 @@
+//! Demosaicing for non-Bayer patterns (Fujifilm X-Trans and others): our own edge-weighted colour-difference
+//! interpolation.
+//!
+//! 1. Green at non-green sites: robust weighted mean of the green samples in the 3×3 neighbourhood (5×5 if
+//!    none), each weighted by `1 / (ε + |g − median|)` so samples across an edge contribute little.
+//! 2. Red and blue everywhere: `G(p) + Σ w_q (C(q) − G(q)) / Σ w_q` over same-colour samples `q` in the 5×5
+//!    neighbourhood, with `w_q = 1 / (d²(p, q) · (ε + |G(q) − G(p)|))` — spatially close samples on the same
+//!    side of an edge dominate.
+
+use super::Mosaic;
+use crate::Rgb32f;
+use lightcraft_raster::par_rows;
+
+const EPS: f32 = 1e-3;
+
+fn green(m: &Mosaic) -> Vec<f32> {
+    let (w, h) = (m.w, m.h);
+    let mut g = vec![0f32; w * h];
+    par_rows(&mut g, w, |y, row| {
+        let y = y as isize;
+        let mut vals: Vec<f32> = Vec::with_capacity(24);
+        for (x, o) in row.iter_mut().enumerate() {
+            let x = x as isize;
+            if m.color(x, y) == 1 {
+                *o = m.get(x, y);
+                continue;
+            }
+            vals.clear();
+            for r in [1isize, 2] {
+                for dy in -r..=r {
+                    for dx in -r..=r {
+                        if (dx != 0 || dy != 0) && m.color(x + dx, y + dy) == 1 {
+                            vals.push(m.get(x + dx, y + dy));
+                        }
+                    }
+                }
+                if !vals.is_empty() {
+                    break;
+                }
+            }
+            if vals.is_empty() {
+                *o = m.get(x, y);
+                continue;
+            }
+            let mut sorted = vals.clone();
+            sorted.sort_by(|a, b| a.total_cmp(b));
+            let med = sorted[sorted.len() / 2];
+            let (mut s, mut ws) = (0.0f32, 0.0f32);
+            for &v in &vals {
+                let wt = 1.0 / (EPS + (v - med).abs());
+                s += wt * v;
+                ws += wt;
+            }
+            *o = s / ws;
+        }
+    });
+    g
+}
+
+pub(crate) fn directional(m: &Mosaic) -> Rgb32f {
+    let (w, h) = (m.w, m.h);
+    let g = green(m);
+    let gg = |x: isize, y: isize| g[super::reflect(y, h) * w + super::reflect(x, w)];
+    let mut out = Rgb32f::new(w, h);
+    par_rows(&mut out.data, w, |y, row| {
+        let y = y as isize;
+        for (x, px) in row.iter_mut().enumerate() {
+            let x = x as isize;
+            let own = m.color(x, y) as usize;
+            let gp = gg(x, y);
+            px[1] = gp;
+            for c in [0usize, 2] {
+                if c == own {
+                    px[c] = m.get(x, y);
+                    continue;
+                }
+                let (mut s, mut ws) = (0.0f32, 0.0f32);
+                for dy in -2isize..=2 {
+                    for dx in -2isize..=2 {
+                        if m.color(x + dx, y + dy) as usize != c {
+                            continue;
+                        }
+                        let gq = gg(x + dx, y + dy);
+                        let d2 = (dx * dx + dy * dy) as f32;
+                        let wt = 1.0 / (d2 * (EPS + (gq - gp).abs()));
+                        s += wt * (m.get(x + dx, y + dy) - gq);
+                        ws += wt;
+                    }
+                }
+                px[c] = if ws > 0.0 { gp + s / ws } else { gp };
+            }
+        }
+    });
+    out
+}
diff --git a/crates/raw/src/dng.rs b/crates/raw/src/dng.rs
new file mode 100644
index 00000000..f25eb604
--- /dev/null
+++ b/crates/raw/src/dng.rs
@@ -0,0 +1,167 @@
+//! DNG (Adobe Digital Negative Specification 1.7): raw IFD selection, pixel data (via [`crate::tiffraw`]),
+//! linearization, black/white levels, active area, default crop, CFA description, colour tags, opcode lists.
+
+use crate::tiffraw::{Packing, read_image};
+use crate::{BlackLevel, Cfa, ColorData, Mat3, RawData, RawError, RawFormat, RawImage, Rect, Result, opcodes};
+use lightcraft_color::Xy;
+use lightcraft_geom::Orientation;
+use lightcraft_tiff::tags::{self as t, photometric};
+use lightcraft_tiff::{Ifd, Tiff};
+
+/// The main raw IFD: full-resolution (NewSubfileType 0) CFA or LinearRaw image with the most pixels.
+pub(crate) fn raw_ifd(tiff: &Tiff) -> Option<&Ifd> {
+    tiff.all_ifds()
+        .into_iter()
+        .filter(|i| i.u32(t::NEW_SUBFILE_TYPE).unwrap_or(0) == 0)
+        .filter(|i| matches!(i.u16(t::PHOTOMETRIC), Some(photometric::CFA) | Some(photometric::LINEAR_RAW)))
+        .max_by_key(|i| i.u64(t::IMAGE_WIDTH).unwrap_or(0) * i.u64(t::IMAGE_LENGTH).unwrap_or(0))
+}
+
+fn mat3(v: Option<Vec<f64>>) -> Option<Mat3> {
+    let v = v?;
+    if v.len() != 9 || v.iter().any(|x| !x.is_finite()) {
+        return None;
+    }
+    let m = Mat3([[v[0], v[1], v[2]], [v[3], v[4], v[5]], [v[6], v[7], v[8]]]);
+    (m.determinant().abs() > 1e-12).then_some(m)
+}
+
+fn vec3(v: Option<Vec<f64>>) -> Option<[f64; 3]> {
+    let v = v?;
+    (v.len() == 3 && v.iter().all(|x| x.is_finite() && *x > 0.0)).then(|| [v[0], v[1], v[2]])
+}
+
+/// Colour tags: DNG puts them in IFD0, but some writers use the raw IFD; prefer the raw IFD.
+pub(crate) fn color_data(ifd0: &Ifd, raw: &Ifd) -> ColorData {
+    let get = |tag: u16| raw.f64s(tag).or_else(|| ifd0.f64s(tag));
+    let geti = |tag: u16| raw.u16(tag).or_else(|| ifd0.u16(tag));
+    ColorData {
+        illuminant: [geti(t::CALIBRATION_ILLUMINANT_1).unwrap_or(0), geti(t::CALIBRATION_ILLUMINANT_2).unwrap_or(0)],
+        color_matrix: [mat3(get(t::COLOR_MATRIX_1)), mat3(get(t::COLOR_MATRIX_2))],
+        forward_matrix: [mat3(get(t::FORWARD_MATRIX_1)), mat3(get(t::FORWARD_MATRIX_2))],
+        camera_calibration: [mat3(get(t::CAMERA_CALIBRATION_1)), mat3(get(t::CAMERA_CALIBRATION_2))],
+        analog_balance: vec3(get(t::ANALOG_BALANCE)),
+        as_shot_neutral: vec3(get(t::AS_SHOT_NEUTRAL)),
+        as_shot_white_xy: get(t::AS_SHOT_WHITE_XY).filter(|v| v.len() == 2 && v[0] > 0.0 && v[1] > 0.0).map(|v| Xy::new(v[0], v[1])),
+        baseline_exposure: get(t::BASELINE_EXPOSURE).and_then(|v| v.first().copied()).filter(|v| v.is_finite()).unwrap_or(0.0)
+            + get(t::BASELINE_EXPOSURE_OFFSET).and_then(|v| v.first().copied()).filter(|v| v.is_finite()).unwrap_or(0.0),
+    }
+}
+
+pub(crate) fn decode(bytes: &[u8]) -> Result<RawImage> {
+    let tiff = Tiff::parse(bytes)?;
+    let ifd0 = &tiff.ifds[0];
+    let raw = raw_ifd(&tiff).ok_or_else(|| RawError::Corrupt("DNG without a raw image IFD".into()))?;
+    let info = raw.image()?;
+    let (w, h, cpp) = (info.width as usize, info.height as usize, info.samples_per_pixel as usize);
+    if !(1..=4).contains(&cpp) {
+        return Err(RawError::Unsupported(format!("{cpp} samples per pixel")));
+    }
+    if info.compression == t::compression::LOSSY_JPEG {
+        return Err(RawError::Unsupported("lossy (DCT) JPEG DNG".into()));
+    }
+    if info.compression == t::compression::JPEG_XL {
+        return Err(RawError::Unsupported("JPEG XL DNG".into()));
+    }
+    let mut data = read_image(bytes, &info, tiff.order, Packing::Msb)?;
+    let bits = info.bits() as u32;
+
+    // linearization table (integer data only)
+    let mut linearized = false;
+    if let (RawData::U16(v), Some(table)) = (&mut data, raw.u64s(t::LINEARIZATION_TABLE).filter(|t| !t.is_empty())) {
+        let table: Vec<u16> = table.iter().map(|&x| x.min(65535) as u16).collect();
+        let last = table.len() - 1;
+        v.iter_mut().for_each(|s| *s = table[(*s as usize).min(last)]);
+        linearized = true;
+    }
+
+    let float = matches!(data, RawData::F32(_));
+    let default_white = if float { 1.0 } else { ((1u64 << bits.min(16)) - 1) as f32 };
+    let mut white: Vec<f32> = raw.f64s(t::WHITE_LEVEL).unwrap_or_default().into_iter().map(|v| v as f32).filter(|v| *v > 0.0).collect();
+    if white.is_empty() {
+        white = vec![default_white];
+    }
+
+    let dim = raw.u64s(t::BLACK_LEVEL_REPEAT_DIM).unwrap_or_default();
+    let (br, bc) = match dim.as_slice() {
+        [r, c] if (1..=16).contains(r) && (1..=16).contains(c) => (*r as usize, *c as usize),
+        _ => (1, 1),
+    };
+    let mut values: Vec<f32> = raw.f64s(t::BLACK_LEVEL).unwrap_or_default().into_iter().map(|v| v as f32).collect();
+    if values.len() != br * bc * cpp {
+        values = vec![values.first().copied().unwrap_or(0.0)];
+    }
+    let (br, bc) = if values.len() == 1 { (1, 1) } else { (br, bc) };
+    let floats = |tag| raw.f64s(tag).unwrap_or_default().into_iter().map(|v| v as f32).collect::<Vec<f32>>();
+    let black =
+        BlackLevel { repeat_rows: br, repeat_cols: bc, values, delta_h: floats(t::BLACK_LEVEL_DELTA_H), delta_v: floats(t::BLACK_LEVEL_DELTA_V) };
+
+    let active_area = match raw.u64s(t::ACTIVE_AREA).as_deref() {
+        Some([top, left, bottom, right]) if bottom > top && right > left && (*bottom as usize) <= h && (*right as usize) <= w => {
+            Rect::new(*left as usize, *top as usize, (*right - *left) as usize, (*bottom - *top) as usize)
+        }
+        _ => Rect::new(0, 0, w, h),
+    };
+    let crop = match (raw.f64s(t::DEFAULT_CROP_ORIGIN).as_deref(), raw.f64s(t::DEFAULT_CROP_SIZE).as_deref()) {
+        (Some([x, y]), Some([cw, ch])) if *cw >= 1.0 && *ch >= 1.0 && *x >= 0.0 && *y >= 0.0 => {
+            Rect::new(x.round() as usize, y.round() as usize, cw.round() as usize, ch.round() as usize).clipped(active_area.width, active_area.height)
+        }
+        _ => Rect::new(0, 0, active_area.width, active_area.height),
+    };
+
+    let cfa = if info.photometric == photometric::CFA {
+        if raw.u16(t::CFA_LAYOUT).unwrap_or(1) != 1 {
+            return Err(RawError::Unsupported("non-rectangular CFA layout".into()));
+        }
+        let dim = raw.u64s(t::CFA_REPEAT_PATTERN_DIM).unwrap_or_else(|| vec![2, 2]);
+        let (rows, cols) = match dim.as_slice() {
+            [r, c] if (1..=16).contains(r) && (1..=16).contains(c) => (*r as usize, *c as usize),
+            _ => return Err(RawError::Corrupt("bad CFARepeatPatternDim".into())),
+        };
+        let pat = raw.bytes(t::CFA_PATTERN_EP).ok_or_else(|| RawError::Corrupt("missing CFAPattern".into()))?;
+        if pat.len() != rows * cols {
+            return Err(RawError::Corrupt("CFAPattern size mismatch".into()));
+        }
+        let planes = raw.bytes(t::CFA_PLANE_COLOR).map(|b| b.to_vec()).unwrap_or_else(|| vec![0, 1, 2]);
+        let pattern: Vec<u8> = pat
+            .iter()
+            .map(|c| planes.iter().position(|p| p == c).map(|i| i as u8).filter(|&i| i < 3))
+            .collect::<Option<_>>()
+            .ok_or_else(|| RawError::Unsupported("CFA with more than three colours".into()))?;
+        if cpp != 1 {
+            return Err(RawError::Unsupported("CFA data with several samples per pixel".into()));
+        }
+        Some(Cfa { width: cols, height: rows, pattern })
+    } else {
+        None
+    };
+
+    let opcode_list = |tag| raw.bytes(tag).map(opcodes::parse_list).unwrap_or_default();
+    let opcodes =
+        crate::OpcodeLists { list1: opcode_list(t::OPCODE_LIST_1), list2: opcode_list(t::OPCODE_LIST_2), list3: opcode_list(t::OPCODE_LIST_3) };
+
+    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::Dng,
+        width: w,
+        height: h,
+        cpp,
+        data,
+        cfa,
+        bits,
+        black,
+        white,
+        active_area,
+        crop,
+        orientation: Orientation::from_exif(ifd0.u16(t::ORIENTATION).unwrap_or(1)),
+        color: color_data(ifd0, raw),
+        wb_multipliers: None,
+        linearized,
+        opcodes,
+        metadata,
+    };
+    img.validate()?;
+    Ok(img)
+}
diff --git a/crates/raw/src/dngwrite.rs b/crates/raw/src/dngwrite.rs
new file mode 100644
index 00000000..6bdeeb27
--- /dev/null
+++ b/crates/raw/src/dngwrite.rs
@@ -0,0 +1,201 @@
+//! Writing DNG files from a [`RawImage`] (uncompressed or lossless-JPEG tiles), used for DNG export and for
+//! bit-exact decoder round-trip tests. The raw image is stored in IFD0 (allowed by DNG 1.7 §"File Structure").
+
+use crate::{RawData, RawError, RawImage, Result, ljpeg, opcodes};
+use lightcraft_color::Mat3;
+use lightcraft_tiff::tags::{self as t, compression, photometric};
+use lightcraft_tiff::writer::{rational, srational};
+use lightcraft_tiff::{ByteOrder, IfdBuilder, ImageData, TiffWriter, Value};
+use rayon::prelude::*;
+
+#[derive(Clone, Copy, Debug, PartialEq, Eq)]
+pub enum DngCompression {
+    Uncompressed,
+    /// Lossless JPEG tiles of `tile × tile` pixels (`tile` even). CFA tiles are encoded as two interleaved
+    /// components of half the tile width (the common DNG layout), linear data with one component per sample.
+    Lj92 {
+        tile: u32,
+    },
+}
+
+#[derive(Clone, Debug, PartialEq)]
+pub struct DngWriteOptions {
+    pub compression: DngCompression,
+    pub order: ByteOrder,
+    /// XMP packet to embed (tag 700); defaults to one generated from the image metadata.
+    pub xmp: Option<String>,
+}
+
+impl Default for DngWriteOptions {
+    fn default() -> Self {
+        Self { compression: DngCompression::Lj92 { tile: 256 }, order: ByteOrder::Little, xmp: None }
+    }
+}
+
+fn srat_mat(m: &Mat3) -> Value {
+    Value::SRational(m.0.iter().flatten().map(|&v| srational(v)).collect())
+}
+
+fn rat_vec(v: &[f64]) -> Value {
+    Value::Rational(v.iter().map(|&x| rational(x)).collect())
+}
+
+/// Serialise `raw` as a DNG.
+pub fn write_dng(raw: &RawImage, opts: &DngWriteOptions) -> Result<Vec<u8>> {
+    raw.validate()?;
+    let (w, h, cpp) = (raw.width, raw.height, raw.cpp);
+    let mut ifd = IfdBuilder::new();
+    ifd.set(t::NEW_SUBFILE_TYPE, Value::Long(vec![0]));
+    ifd.set(t::IMAGE_WIDTH, Value::Long(vec![w as u32]));
+    ifd.set(t::IMAGE_LENGTH, Value::Long(vec![h as u32]));
+    ifd.set(t::SAMPLES_PER_PIXEL, Value::Short(vec![cpp as u16]));
+    ifd.set(t::PLANAR_CONFIGURATION, Value::Short(vec![1]));
+    ifd.set(t::PHOTOMETRIC, Value::Short(vec![if raw.cfa.is_some() { photometric::CFA } else { photometric::LINEAR_RAW }]));
+    ifd.set(t::DNG_VERSION, Value::Byte(vec![1, 4, 0, 0]));
+    ifd.set(t::DNG_BACKWARD_VERSION, Value::Byte(vec![1, 3, 0, 0]));
+    let make = raw.metadata.make.clone().unwrap_or_default();
+    let model = raw.metadata.model.clone().unwrap_or_default();
+    if !make.is_empty() {
+        ifd.set(t::MAKE, Value::Ascii(make.clone()));
+    }
+    if !model.is_empty() {
+        ifd.set(t::MODEL, Value::Ascii(model.clone()));
+    }
+    let unique = format!("{make} {model}").trim().to_string();
+    ifd.set(t::UNIQUE_CAMERA_MODEL, Value::Ascii(if unique.is_empty() { "LightCraft".into() } else { unique }));
+    ifd.set(t::ORIENTATION, Value::Short(vec![raw.orientation.to_exif()]));
+    ifd.set(t::SOFTWARE, Value::Ascii("LightCraft".into()));
+    let xmp = opts.xmp.clone().unwrap_or_else(|| lightcraft_meta::write_xmp(&raw.metadata, None));
+    ifd.set(t::XMP, Value::Byte(xmp.into_bytes()));
+
+    if let Some(cfa) = &raw.cfa {
+        ifd.set(t::CFA_REPEAT_PATTERN_DIM, Value::Short(vec![cfa.height as u16, cfa.width as u16]));
+        ifd.set(t::CFA_PATTERN_EP, Value::Byte(cfa.pattern.clone()));
+        ifd.set(t::CFA_PLANE_COLOR, Value::Byte(vec![0, 1, 2]));
+        ifd.set(t::CFA_LAYOUT, Value::Short(vec![1]));
+    }
+    let b = &raw.black;
+    ifd.set(t::BLACK_LEVEL_REPEAT_DIM, Value::Short(vec![b.repeat_rows as u16, b.repeat_cols as u16]));
+    let mut bl: Vec<f64> = b.values.iter().map(|&v| v as f64).collect();
+    if bl.len() == 1 && cpp > 1 {
+        bl = vec![bl[0]; cpp];
+    }
+    ifd.set(t::BLACK_LEVEL, rat_vec(&bl));
+    if !b.delta_h.is_empty() {
+        ifd.set(t::BLACK_LEVEL_DELTA_H, Value::SRational(b.delta_h.iter().map(|&v| srational(v as f64)).collect()));
+    }
+    if !b.delta_v.is_empty() {
+        ifd.set(t::BLACK_LEVEL_DELTA_V, Value::SRational(b.delta_v.iter().map(|&v| srational(v as f64)).collect()));
+    }
+    let white: Vec<u32> = (0..cpp).map(|s| raw.white_at(s).round().clamp(0.0, u32::MAX as f32) as u32).collect();
+    if matches!(raw.data, RawData::F32(_)) {
+        ifd.set(t::WHITE_LEVEL, rat_vec(&(0..cpp).map(|s| raw.white_at(s) as f64).collect::<Vec<_>>()));
+    } else {
+        ifd.set(t::WHITE_LEVEL, Value::Long(white));
+    }
+    let a = raw.active_area;
+    ifd.set(t::ACTIVE_AREA, Value::Long(vec![a.y as u32, a.x as u32, (a.y + a.height) as u32, (a.x + a.width) as u32]));
+    ifd.set(t::DEFAULT_CROP_ORIGIN, Value::Long(vec![raw.crop.x as u32, raw.crop.y as u32]));
+    ifd.set(t::DEFAULT_CROP_SIZE, Value::Long(vec![raw.crop.width as u32, raw.crop.height as u32]));
+
+    let c = &raw.color;
+    for (i, tag) in [t::CALIBRATION_ILLUMINANT_1, t::CALIBRATION_ILLUMINANT_2].into_iter().enumerate() {
+        if c.color_matrix[i].is_some() {
+            ifd.set(tag, Value::Short(vec![c.illuminant[i]]));
+        }
+    }
+    for (pair, tags) in [
+        (&c.color_matrix, [t::COLOR_MATRIX_1, t::COLOR_MATRIX_2]),
+        (&c.forward_matrix, [t::FORWARD_MATRIX_1, t::FORWARD_MATRIX_2]),
+        (&c.camera_calibration, [t::CAMERA_CALIBRATION_1, t::CAMERA_CALIBRATION_2]),
+    ] {
+        for i in 0..2 {
+            if let Some(m) = &pair[i] {
+                ifd.set(tags[i], srat_mat(m));
+            }
+        }
+    }
+    if let Some(v) = c.analog_balance {
+        ifd.set(t::ANALOG_BALANCE, rat_vec(&v));
+    }
+    if let Some(v) = c.as_shot_neutral {
+        ifd.set(t::AS_SHOT_NEUTRAL, rat_vec(&v));
+    } else if let Some(xy) = c.as_shot_white_xy {
+        ifd.set(t::AS_SHOT_WHITE_XY, rat_vec(&[xy.x, xy.y]));
+    }
+    ifd.set(t::BASELINE_EXPOSURE, Value::SRational(vec![srational(c.baseline_exposure)]));
+    for (list, tag) in [(&raw.opcodes.list1, t::OPCODE_LIST_1), (&raw.opcodes.list2, t::OPCODE_LIST_2), (&raw.opcodes.list3, t::OPCODE_LIST_3)] {
+        if !list.is_empty() {
+            ifd.set(tag, Value::Undefined(opcodes::write_list(list)));
+        }
+    }
+
+    match &raw.data {
+        RawData::F32(v) => {
+            ifd.set(t::BITS_PER_SAMPLE, Value::Short(vec![32; cpp]));
+            ifd.set(t::SAMPLE_FORMAT, Value::Short(vec![3; cpp]));
+            ifd.set(t::COMPRESSION, Value::Short(vec![compression::NONE]));
+            let mut bytes = Vec::with_capacity(v.len() * 4);
+            for &x in v {
+                let mut tmp = Vec::with_capacity(4);
+                opts.order.put_u32(&mut tmp, x.to_bits());
+                bytes.extend_from_slice(&tmp);
+            }
+            ifd.set_image(ImageData::Strips { rows_per_strip: h as u32, strips: vec![bytes] });
+        }
+        RawData::U16(v) => {
+            let maxv = v.iter().copied().max().unwrap_or(0);
+            let bits = (16 - maxv.leading_zeros()).max(raw.bits.min(16)).max(2) as u8;
+            ifd.set(
+                t::BITS_PER_SAMPLE,
+                Value::Short(vec![if matches!(opts.compression, DngCompression::Uncompressed) { 16 } else { bits as u16 }; cpp]),
+            );
+            match opts.compression {
+                DngCompression::Uncompressed => {
+                    ifd.set(t::COMPRESSION, Value::Short(vec![compression::NONE]));
+                    let rows = (65536 / (w * cpp * 2)).clamp(1, h);
+                    let strips = v
+                        .chunks(rows * w * cpp)
+                        .map(|s| {
+                            let mut b = Vec::with_capacity(s.len() * 2);
+                            s.iter().for_each(|&x| opts.order.put_u16(&mut b, x));
+                            b
+                        })
+                        .collect();
+                    ifd.set_image(ImageData::Strips { rows_per_strip: rows as u32, strips });
+                }
+                DngCompression::Lj92 { tile } => {
+                    let tile = (tile.max(2) & !1) as usize;
+                    if tile > 65534 {
+                        return Err(RawError::Limit("tile too large"));
+                    }
+                    ifd.set(t::COMPRESSION, Value::Short(vec![compression::JPEG]));
+                    let (ta, td) = (w.div_ceil(tile), h.div_ceil(tile));
+                    let tiles: Vec<Vec<u8>> = (0..ta * td)
+                        .into_par_iter()
+                        .map(|i| {
+                            let (tx, ty) = (i % ta * tile, i / ta * tile);
+                            let mut buf = Vec::with_capacity(tile * tile * cpp);
+                            for y in 0..tile {
+                                let sy = (ty + y).min(h - 1);
+                                for x in 0..tile {
+                                    let sx = (tx + x).min(w - 1);
+                                    for s in 0..cpp {
+                                        buf.push(v[(sy * w + sx) * cpp + s]);
+                                    }
+                                }
+                            }
+                            if cpp == 1 {
+                                ljpeg::encode(&buf, tile / 2, tile, 2, bits, 1, 0)
+                            } else {
+                                ljpeg::encode(&buf, tile, tile, cpp, bits, 1, 0)
+                            }
+                        })
+                        .collect();
+                    ifd.set_image(ImageData::Tiles { tile_width: tile as u32, tile_height: tile as u32, tiles });
+                }
+            }
+        }
+    }
+    Ok(TiffWriter::new(opts.order, false).write(&[ifd])?)
+}
diff --git a/crates/raw/src/highlight.rs b/crates/raw/src/highlight.rs
new file mode 100644
index 00000000..388174ee
--- /dev/null
+++ b/crates/raw/src/highlight.rs
@@ -0,0 +1,180 @@
+//! Clip-aware highlight handling on demosaiced camera RGB (before white balance, white level = 1.0).
+//!
+//! - [`clip_neutral`]: white-balance-aware clipping so that sensor-clipped areas render neutral instead of
+//!   magenta/cyan (each channel limited to the smallest white-balanced clip level).
+//! - [`reconstruct`]: where only some channels are clipped, rebuild them from the unclipped channels using the
+//!   local chromaticity of nearby unclipped pixels (diffused into the clipped region with a coarse-to-fine
+//!   normalised-convolution fill). Fully clipped pixels become neutral at the brightest plausible level.
+
+use lightcraft_raster::Rgb32f;
+
+/// Clip every channel of `wb ⊙ img` at `min_c(wb_c) · clip` — i.e. at the lowest channel's clip level after WB —
+/// then divide the multipliers back out. `clip` is the sensor clip in normalised units (≈ 1.0).
+pub fn clip_neutral(img: &mut Rgb32f, wb: [f32; 3], clip: f32) {
+    let limit = wb.iter().cloned().fold(f32::MAX, f32::min) * clip;
+    for p in &mut img.data {
+        for c in 0..3 {
+            p[c] = (p[c] * wb[c]).min(limit) / wb[c];
+        }
+    }
+}
+
+/// Fill `values` (per-pixel vectors) where `valid` is false from valid neighbours, coarse to fine.
+fn fill_invalid(w: usize, h: usize, values: &mut [[f32; 3]], valid: &mut [bool]) {
+    if w == 0 || h == 0 || valid.iter().all(|&v| v) || !valid.iter().any(|&v| v) {
+        return;
+    }
+    // Build a pyramid of weighted means, fill each level from the next coarser one.
+    let (cw, ch) = (w.div_ceil(2), h.div_ceil(2));
+    let mut cv = vec![[0f32; 3]; cw * ch];
+    let mut cval = vec![false; cw * ch];
+    for y in 0..ch {
+        for x in 0..cw {
+            let (mut s, mut n) = ([0f32; 3], 0);
+            for dy in 0..2 {
+                for dx in 0..2 {
+                    let (sx, sy) = (2 * x + dx, 2 * y + dy);
+                    if sx < w && sy < h && valid[sy * w + sx] {
+                        for c in 0..3 {
+                            s[c] += values[sy * w + sx][c];
+                        }
+                        n += 1;
+                    }
+                }
+            }
+            if n > 0 {
+                cv[y * cw + x] = s.map(|v| v / n as f32);
+                cval[y * cw + x] = true;
+            }
+        }
+    }
+    if cw * ch < w * h {
+        fill_invalid(cw, ch, &mut cv, &mut cval);
+    } else {
+        // cannot shrink further (1×1): nothing valid anywhere handled above
+        return;
+    }
+    for y in 0..h {
+        for x in 0..w {
+            let i = y * w + x;
+            if !valid[i] {
+                // bilinear upsample of the coarse level
+                let fx = ((x as f32 + 0.5) / 2.0 - 0.5).clamp(0.0, (cw - 1) as f32);
+                let fy = ((y as f32 + 0.5) / 2.0 - 0.5).clamp(0.0, (ch - 1) as f32);
+                let (x0, y0) = (fx.floor() as usize, fy.floor() as usize);
+                let (x1, y1) = ((x0 + 1).min(cw - 1), (y0 + 1).min(ch - 1));
+                let (tx, ty) = (fx - x0 as f32, fy - y0 as f32);
+                let g = |xx: usize, yy: usize| cv[yy * cw + xx];
+                let mut v = [0f32; 3];
+                for c in 0..3 {
+                    let top = g(x0, y0)[c] + (g(x1, y0)[c] - g(x0, y0)[c]) * tx;
+                    let bot = g(x0, y1)[c] + (g(x1, y1)[c] - g(x0, y1)[c]) * tx;
+                    v[c] = top + (bot - top) * ty;
+                }
+                values[i] = v;
+                valid[i] = true;
+            }
+        }
+    }
+}
+
+/// Reconstruct partially clipped channels. `clip` is the sensor clip level in normalised units (use slightly
+/// below 1.0, e.g. 0.99), `wb` the white-balance multipliers that will be applied afterwards.
+/// Returns the number of pixels that had at least one clipped channel.
+pub fn reconstruct(img: &mut Rgb32f, wb: [f32; 3], clip: f32) -> usize {
+    let (w, h) = (img.width, img.height);
+    let n = w * h;
+    let clipped: Vec<[bool; 3]> = img.data.iter().map(|p| [p[0] >= clip, p[1] >= clip, p[2] >= clip]).collect();
+    let count = clipped.iter().filter(|c| c.iter().any(|&b| b)).count();
+    if count == 0 {
+        return 0;
+    }
+    // chromaticity (white-balanced ratios to the channel mean) of unclipped pixels
+    let mut chroma = vec![[1f32 / 3.0; 3]; n];
+    let mut valid = vec![false; n];
+    for i in 0..n {
+        if clipped[i].iter().any(|&b| b) {
+            continue;
+        }
+        let p = img.data[i];
+        let q = [p[0] * wb[0], p[1] * wb[1], p[2] * wb[2]];
+        let s = q[0] + q[1] + q[2];
+        if s > 1e-4 {
+            chroma[i] = q.map(|v| v.max(0.0) / s);
+            valid[i] = true;
+        }
+    }
+    fill_invalid(w, h, &mut chroma, &mut valid);
+    let max_level = wb.iter().cloned().fold(0.0f32, f32::max) * clip;
+    for i in 0..n {
+        let cl = clipped[i];
+        if !cl.iter().any(|&b| b) {
+            continue;
+        }
+        let p = img.data[i];
+        let q = [p[0] * wb[0], p[1] * wb[1], p[2] * wb[2]];
+        let r = chroma[i];
+        // estimate the white-balanced channel sum from the unclipped channels
+        let ests: Vec<f32> = (0..3).filter(|&c| !cl[c] && r[c] > 1e-3).map(|c| q[c] / r[c]).collect();
+        let mut out = q;
+        if ests.is_empty() {
+            let v = q.iter().cloned().fold(0.0f32, f32::max).max(max_level);
+            out = [v; 3];
+        } else {
+            let sum = ests.iter().sum::<f32>() / ests.len() as f32;
+            for c in 0..3 {
+                if cl[c] {
+                    out[c] = q[c].max(sum * r[c]);
+                }
+            }
+        }
+        img.data[i] = [out[0] / wb[0], out[1] / wb[1], out[2] / wb[2]];
+    }
+    count
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn clip_neutral_makes_clipped_white_neutral() {
+        let wb = [2.0, 1.0, 1.5];
+        let mut img = Rgb32f::filled(2, 1, [1.0, 1.0, 1.0]);
+        img.data[1] = [0.2, 0.3, 0.4];
+        clip_neutral(&mut img, wb, 1.0);
+        let p = img.data[0];
+        let q = [p[0] * wb[0], p[1] * wb[1], p[2] * wb[2]];
+        assert!((q[0] - q[1]).abs() < 1e-6 && (q[1] - q[2]).abs() < 1e-6);
+        assert_eq!(img.data[1], [0.2, 0.3, 0.4]);
+    }
+
+    #[test]
+    fn reconstruct_restores_clipped_channel() {
+        // a warm gradient whose red channel clips in the bright half
+        let truth = Rgb32f::from_fn(64, 16, |x, _| {
+            let v = 0.3 + 1.2 * x as f32 / 63.0;
+            [v * 1.0, v * 0.6, v * 0.35]
+        });
+        let mut img = truth.map(|p| p.map(|v| v.min(1.0)));
+        let before: f32 = img.data.iter().zip(&truth.data).map(|(a, b)| (a[0] - b[0]).abs()).sum();
+        let n = reconstruct(&mut img, [1.0, 1.0, 1.0], 0.999);
+        assert!(n > 0);
+        let after: f32 = img.data.iter().zip(&truth.data).map(|(a, b)| (a[0] - b[0]).abs()).sum();
+        assert!(after < before * 0.2, "before {before} after {after}");
+        // unclipped pixels untouched
+        assert_eq!(img.get(0, 0), truth.get(0, 0));
+    }
+
+    #[test]
+    fn fully_clipped_and_no_clipping() {
+        let mut img = Rgb32f::filled(4, 4, [1.0; 3]);
+        assert_eq!(reconstruct(&mut img, [2.0, 1.0, 1.5], 0.99), 16);
+        let p = img.data[0];
+        assert!((p[0] * 2.0 - p[1]).abs() < 1e-5 && (p[2] * 1.5 - p[1]).abs() < 1e-5);
+        let mut img = Rgb32f::filled(4, 4, [0.5; 3]);
+        assert_eq!(reconstruct(&mut img, [2.0, 1.0, 1.5], 0.99), 0);
+        let mut empty = Rgb32f::new(0, 0);
+        assert_eq!(reconstruct(&mut empty, [1.0; 3], 0.99), 0);
+    }
+}
diff --git a/crates/raw/src/lib.rs b/crates/raw/src/lib.rs
index 45278f22..cde39418 100644
--- a/crates/raw/src/lib.rs
+++ b/crates/raw/src/lib.rs
@@ -1 +1,490 @@
+//! RAW decoding for LightCraft.
+//!
+//! - [`probe`] recognises raw containers; [`decode`] turns a file into a [`RawImage`] (sensor data + everything
+//!   needed to render it: CFA, black/white levels, active area, default crop, orientation, DNG colour tags,
+//!   opcode lists, [`Metadata`]); [`embedded_preview`] returns the largest embedded JPEG.
+//! - [`RawImage::normalized`] subtracts black, scales white to 1.0 and crops to the active area (applying DNG
+//!   `OpcodeList1`/`OpcodeList2`); [`demosaic`] turns CFA data into camera-RGB [`Rgb32f`];
+//!   [`RawImage::develop`] does all of it plus `OpcodeList3` and the default crop.
+//! - [`color`] implements the DNG colour model (dual-illuminant interpolation, forward matrices, white balance)
+//!   and produces camera → linear Rec.2020 D65 matrices.
+//!
+//! Formats: DNG (uncompressed, lossless JPEG, Deflate incl. floating point, tiled/stripped, CFA and LinearRaw),
+//! Canon CR2, Nikon NEF (uncompressed), Sony ARW (uncompressed, lossless). See the crate README for sources and
+//! gaps. The decoders never panic on malformed input.
 #![forbid(unsafe_code)]
+
+pub mod color;
+pub mod demosaic;
+mod dng;
+pub mod dngwrite;
+pub mod highlight;
+pub mod ljpeg;
+pub mod opcodes;
+mod preview;
+mod tiffraw;
+mod unpack;
+mod vendor;
+
+pub use demosaic::{Method, demosaic};
+pub use dngwrite::{DngCompression, DngWriteOptions, write_dng};
+pub use lightcraft_color::Mat3;
+pub use lightcraft_geom::Orientation;
+pub use lightcraft_meta::Metadata;
+pub use lightcraft_raster::Rgb32f;
+pub use opcodes::{Opcode, OpcodeLists};
+pub use preview::embedded_preview;
+
+use lightcraft_color::Xy;
+use lightcraft_tiff::{Tiff, TiffError};
+use serde::{Deserialize, Serialize};
+
+/// Upper bound on decoded samples (guards allocations driven by header values).
+pub const MAX_SAMPLES: usize = 1 << 30;
+
+/// Errors from raw decoding.
+#[derive(Debug, Clone, PartialEq, Eq, thiserror::Error)]
+pub enum RawError {
+    #[error("not a recognised raw file")]
+    NotRaw,
+    #[error("unsupported raw variant: {0}")]
+    Unsupported(String),
+    #[error("corrupt raw data: {0}")]
+    Corrupt(String),
+    #[error("limit exceeded: {0}")]
+    Limit(&'static str),
+    #[error(transparent)]
+    Tiff(#[from] TiffError),
+}
+
+pub type Result<T> = std::result::Result<T, RawError>;
+
+/// Raw container formats recognised by [`probe`].
+#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash, Serialize, Deserialize)]
+pub enum RawFormat {
+    Dng,
+    Cr2,
+    Cr3,
+    Nef,
+    Nrw,
+    Arw,
+    Raf,
+    Orf,
+    Rw2,
+    Pef,
+    Srw,
+    /// Another TIFF-based raw (3FR, IIQ, ERF, KDC, DCR, MOS, …).
+    OtherTiff,
+}
+
+impl RawFormat {
+    /// Whether [`decode`] supports this container (possibly not every compression inside it).
+    pub fn is_supported(self) -> bool {
+        matches!(self, RawFormat::Dng | RawFormat::Cr2 | RawFormat::Nef | RawFormat::Nrw | RawFormat::Arw)
+    }
+}
+
+/// Recognise a raw file from its first bytes / TIFF structure.
+pub fn probe(bytes: &[u8]) -> Option<RawFormat> {
+    if bytes.len() >= 12 && &bytes[4..12] == b"ftypcrx " {
+        return Some(RawFormat::Cr3);
+    }
+    if bytes.starts_with(b"FUJIFILMCCD-RAW") {
+        return Some(RawFormat::Raf);
+    }
+    if bytes.starts_with(b"IIRO") || bytes.starts_with(b"IIRS") || bytes.starts_with(b"MMOR") {
+        return Some(RawFormat::Orf);
+    }
+    if bytes.starts_with(b"IIU\0") {
+        return Some(RawFormat::Rw2);
+    }
+    if bytes.len() >= 10 && bytes.starts_with(b"II*\0") && &bytes[8..10] == b"CR" {
+        return Some(RawFormat::Cr2);
+    }
+    let (_, _) = Tiff::sniff(bytes)?;
+    let opts = lightcraft_tiff::ParseOptions { max_ifds: 256, ..Default::default() };
+    let t = Tiff::parse_with(bytes, &opts).ok()?;
+    let ifd0 = t.ifds.first()?;
+    if ifd0.contains(lightcraft_tiff::tags::DNG_VERSION) {
+        return Some(RawFormat::Dng);
+    }
+    let make = t.find(lightcraft_tiff::tags::MAKE).and_then(|e| e.value.as_str()).unwrap_or_default().to_ascii_uppercase();
+    let has_cfa = t.all_ifds().iter().any(|i| {
+        i.u16(lightcraft_tiff::tags::PHOTOMETRIC) == Some(lightcraft_tiff::tags::photometric::CFA)
+            || i.u16(lightcraft_tiff::tags::COMPRESSION).is_some_and(|c| c == 34713 || c == 32767 || c == 32769 || c == 32770)
+    });
+    if make.starts_with("NIKON") {
+        return Some(if has_cfa || t.all_ifds().len() > 1 { RawFormat::Nef } else { RawFormat::Nrw });
+    }
+    if make.starts_with("SONY") {
+        return Some(RawFormat::Arw);
+    }
+    if make.starts_with("PENTAX") || make.starts_with("RICOH") {
+        return Some(RawFormat::Pef);
+    }
+    if make.starts_with("SAMSUNG") {
+        return Some(RawFormat::Srw);
+    }
+    if has_cfa {
+        return Some(RawFormat::OtherTiff);
+    }
+    None
+}
+
+/// Decode a raw file.
+pub fn decode(bytes: &[u8]) -> Result<RawImage> {
+    match probe(bytes).ok_or(RawError::NotRaw)? {
+        RawFormat::Dng => dng::decode(bytes),
+        RawFormat::Cr2 => vendor::cr2::decode(bytes),
+        RawFormat::Nef | RawFormat::Nrw => vendor::nef::decode(bytes),
+        RawFormat::Arw => vendor::arw::decode(bytes),
+        other => Err(RawError::Unsupported(format!("{other:?} files are not decoded yet"))),
+    }
+}
+
+/// Sensor samples.
+#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
+pub enum RawData {
+    U16(Vec<u16>),
+    /// Floating-point DNGs (already linear, typically white = 1.0).
+    F32(Vec<f32>),
+}
+
+impl RawData {
+    pub fn len(&self) -> usize {
+        match self {
+            RawData::U16(v) => v.len(),
+            RawData::F32(v) => v.len(),
+        }
+    }
+    pub fn is_empty(&self) -> bool {
+        self.len() == 0
+    }
+    #[inline]
+    pub fn get(&self, i: usize) -> f32 {
+        match self {
+            RawData::U16(v) => v[i] as f32,
+            RawData::F32(v) => v[i],
+        }
+    }
+}
+
+/// Colour filter array pattern. Colour indices: 0 = red, 1 = green, 2 = blue.
+#[derive(Clone, Debug, PartialEq, Eq, Hash, Serialize, Deserialize)]
+pub struct Cfa {
+    pub width: usize,
+    pub height: usize,
+    /// `width × height` colour indices, row-major, anchored at image pixel (0, 0).
+    pub pattern: Vec<u8>,
+}
+
+impl Cfa {
+    /// A 2×2 Bayer pattern from a string such as `"RGGB"`.
+    pub fn bayer(s: &str) -> Option<Cfa> {
+        let p: Vec<u8> = s
+            .chars()
+            .map(|c| {
+                Some(match c {
+                    'R' => 0,
+                    'G' => 1,
+                    'B' => 2,
+                    _ => return None,
+                })
+            })
+            .collect::<Option<_>>()?;
+        (p.len() == 4).then_some(Cfa { width: 2, height: 2, pattern: p })
+    }
+    /// The Fujifilm X-Trans 6×6 layout (as commonly documented), anchored at (0, 0).
+    pub fn xtrans() -> Cfa {
+        let rows = ["GGRGGB", "GGBGGR", "BRGRBG", "GGBGGR", "GGRGGB", "RBGBRG"];
+        let pattern = rows
+            .iter()
+            .flat_map(|r| {
+                r.chars().map(|c| match c {
+                    'R' => 0,
+                    'G' => 1,
+                    _ => 2,
+                })
+            })
+            .collect();
+        Cfa { width: 6, height: 6, pattern }
+    }
+    #[inline]
+    pub fn color_at(&self, x: usize, y: usize) -> u8 {
+        self.pattern[(y % self.height) * self.width + (x % self.width)]
+    }
+    /// The same pattern anchored at `(dx, dy)` of the current anchor (e.g. after cropping to an active area).
+    pub fn shifted(&self, dx: usize, dy: usize) -> Cfa {
+        let pattern = (0..self.height).flat_map(|y| (0..self.width).map(move |x| (x, y))).map(|(x, y)| self.color_at(x + dx, y + dy)).collect();
+        Cfa { width: self.width, height: self.height, pattern }
+    }
+    pub fn is_bayer(&self) -> bool {
+        if self.width != 2 || self.height != 2 {
+            return false;
+        }
+        let mut c = [0; 3];
+        for &p in &self.pattern {
+            if p > 2 {
+                return false;
+            }
+            c[p as usize] += 1;
+        }
+        c == [1, 2, 1]
+    }
+    /// Pattern name like `"RGGB"` for 2×2 patterns.
+    pub fn name(&self) -> String {
+        self.pattern.iter().map(|&c| ['R', 'G', 'B', '?'][c.min(3) as usize]).collect()
+    }
+    fn valid(&self) -> bool {
+        self.width > 0 && self.height > 0 && self.width <= 16 && self.height <= 16 && self.pattern.len() == self.width * self.height
+    }
+}
+
+/// Integer rectangle (pixels).
+#[derive(Clone, Copy, Debug, Default, PartialEq, Eq, Hash, Serialize, Deserialize)]
+pub struct Rect {
+    pub x: usize,
+    pub y: usize,
+    pub width: usize,
+    pub height: usize,
+}
+
+impl Rect {
+    pub fn new(x: usize, y: usize, width: usize, height: usize) -> Rect {
+        Rect { x, y, width, height }
+    }
+    /// Clip to a `w × h` area.
+    pub fn clipped(self, w: usize, h: usize) -> Rect {
+        let x = self.x.min(w);
+        let y = self.y.min(h);
+        Rect { x, y, width: self.width.min(w - x), height: self.height.min(h - y) }
+    }
+}
+
+/// Black level model (DNG `BlackLevelRepeatDim`, `BlackLevel`, `BlackLevelDeltaH/V`). The repeat pattern and
+/// the deltas are anchored at the active area's top-left corner.
+#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
+pub struct BlackLevel {
+    pub repeat_rows: usize,
+    pub repeat_cols: usize,
+    /// `repeat_rows × repeat_cols × cpp` values.
+    pub values: Vec<f32>,
+    /// Per active-area column.
+    pub delta_h: Vec<f32>,
+    /// Per active-area row.
+    pub delta_v: Vec<f32>,
+}
+
+impl Default for BlackLevel {
+    fn default() -> Self {
+        BlackLevel { repeat_rows: 1, repeat_cols: 1, values: vec![0.0], delta_h: vec![], delta_v: vec![] }
+    }
+}
+
+impl BlackLevel {
+    pub fn uniform(v: f32) -> BlackLevel {
+        BlackLevel { values: vec![v], ..Default::default() }
+    }
+    /// Black level at active-area coordinates (x, y) for sample `s` of `cpp`.
+    #[inline]
+    pub fn at(&self, x: usize, y: usize, s: usize, cpp: usize) -> f32 {
+        let (r, c) = (self.repeat_rows.max(1), self.repeat_cols.max(1));
+        let idx = ((y % r) * c + (x % c)) * cpp + s;
+        let base = self.values.get(idx).or_else(|| self.values.first()).copied().unwrap_or(0.0);
+        base + self.delta_h.get(x).copied().unwrap_or(0.0) + self.delta_v.get(y).copied().unwrap_or(0.0)
+    }
+    /// Mean of the repeat-pattern values (for display / heuristics).
+    pub fn mean(&self) -> f32 {
+        if self.values.is_empty() { 0.0 } else { self.values.iter().sum::<f32>() / self.values.len() as f32 }
+    }
+}
+
+/// DNG colour tags (DNG 1.7 chapter 6). Matrices are 3×3 (three-colour cameras).
+#[derive(Clone, Debug, Default, PartialEq, Serialize, Deserialize)]
+pub struct ColorData {
+    /// `CalibrationIlluminant1/2` (Exif LightSource codes; 0 = unknown).
+    pub illuminant: [u16; 2],
+    /// `ColorMatrix1/2`: XYZ → reference camera.
+    pub color_matrix: [Option<Mat3>; 2],
+    /// `ForwardMatrix1/2`: white-balanced camera → XYZ D50.
+    pub forward_matrix: [Option<Mat3>; 2],
+    /// `CameraCalibration1/2`.
+    pub camera_calibration: [Option<Mat3>; 2],
+    pub analog_balance: Option<[f64; 3]>,
+    pub as_shot_neutral: Option<[f64; 3]>,
+    pub as_shot_white_xy: Option<Xy>,
+    /// EV to add for a "normal" rendering (`BaselineExposure` + `BaselineExposureOffset`).
+    pub baseline_exposure: f64,
+}
+
+/// A decoded raw image.
+#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
+pub struct RawImage {
+    pub format: RawFormat,
+    /// Full sensor data dimensions (including masked borders).
+    pub width: usize,
+    pub height: usize,
+    /// Samples per pixel: 1 = CFA or monochrome, 3 = linear RGB (LinearRaw / demosaiced DNG).
+    pub cpp: usize,
+    pub data: RawData,
+    /// CFA pattern anchored at pixel (0, 0) of the full data; `None` for linear / monochrome data.
+    pub cfa: Option<Cfa>,
+    /// Significant bits of the stored samples (informational).
+    pub bits: u32,
+    pub black: BlackLevel,
+    /// White (clip) level per sample (after linearization).
+    pub white: Vec<f32>,
+    /// Area of valid image data within the full data.
+    pub active_area: Rect,
+    /// Default crop relative to the active area.
+    pub crop: Rect,
+    pub orientation: Orientation,
+    pub color: ColorData,
+    /// As-shot white-balance multipliers from vendor maker notes (green = 1), when known (non-DNG).
+    pub wb_multipliers: Option<[f32; 3]>,
+    /// Whether a `LinearizationTable` was applied to `data`.
+    pub linearized: bool,
+    pub opcodes: OpcodeLists,
+    pub metadata: Metadata,
+}
+
+/// Black-subtracted, white-normalised data covering the active area.
+#[derive(Clone, Debug, PartialEq)]
+pub struct Normalized {
+    pub width: usize,
+    pub height: usize,
+    pub cpp: usize,
+    /// Nominal range [0, 1] (values may exceed 1 above white or be negative from noise).
+    pub data: Vec<f32>,
+    /// CFA anchored at (0, 0) of this buffer.
+    pub cfa: Option<Cfa>,
+}
+
+impl RawImage {
+    /// Validate internal consistency (dimensions vs data length, CFA shape, rectangles).
+    pub fn validate(&self) -> Result<()> {
+        if self.width == 0 || self.height == 0 || !(1..=4).contains(&self.cpp) {
+            return Err(RawError::Corrupt("bad dimensions".into()));
+        }
+        if self.data.len() != self.width * self.height * self.cpp {
+            return Err(RawError::Corrupt("data length mismatch".into()));
+        }
+        if let Some(c) = &self.cfa
+            && (!c.valid() || c.pattern.iter().any(|&p| p > 2))
+        {
+            return Err(RawError::Unsupported(format!("CFA pattern {}", c.name())));
+        }
+        let a = self.active_area;
+        if a.width == 0 || a.height == 0 || a.x + a.width > self.width || a.y + a.height > self.height {
+            return Err(RawError::Corrupt("active area outside image".into()));
+        }
+        Ok(())
+    }
+
+    /// White level for sample `s`.
+    pub fn white_at(&self, s: usize) -> f32 {
+        self.white.get(s).or_else(|| self.white.first()).copied().unwrap_or(65535.0)
+    }
+
+    /// Black-subtract, white-scale and crop to the active area, applying `OpcodeList1` (on raw values) first and
+    /// `OpcodeList2` (on normalised values) after. 1.0 = white level.
+    pub fn normalized(&self) -> Result<Normalized> {
+        self.validate()?;
+        let (w, cpp) = (self.width, self.cpp);
+        let a = self.active_area;
+        // stage 1 (full image) opcodes need a mutable float copy only when present
+        let stage1: Option<Vec<f32>> = if self.opcodes.list1.is_empty() {
+            None
+        } else {
+            let mut v: Vec<f32> = (0..self.data.len()).map(|i| self.data.get(i)).collect();
+            opcodes::apply_list(&self.opcodes.list1, &mut v, self.width, self.height, cpp, self.cfa.as_ref(), 65535.0);
+            Some(v)
+        };
+        let src = |i: usize| match &stage1 {
+            Some(v) => v[i],
+            None => self.data.get(i),
+        };
+        let scale: Vec<f32> = (0..cpp)
+            .map(|s| {
+                let range = self.white_at(s) - self.black.mean();
+                if range > 0.0 { 1.0 / range } else { 1.0 }
+            })
+            .collect();
+        let mut out = vec![0f32; a.width * a.height * cpp];
+        let row_len = a.width * cpp;
+        use rayon::prelude::*;
+        out.par_chunks_mut(row_len).enumerate().for_each(|(y, row)| {
+            let base = ((a.y + y) * w + a.x) * cpp;
+            for x in 0..a.width {
+                for s in 0..cpp {
+                    let v = src(base + x * cpp + s);
+                    let b = self.black.at(x, y, s, cpp);
+                    row[x * cpp + s] = (v - b) * scale[s];
+                }
+            }
+        });
+        let cfa = self.cfa.as_ref().map(|c| c.shifted(a.x, a.y));
+        opcodes::apply_list(&self.opcodes.list2, &mut out, a.width, a.height, cpp, cfa.as_ref(), 1.0);
+        Ok(Normalized { width: a.width, height: a.height, cpp, data: out, cfa })
+    }
+
+    /// Full "raw → camera RGB" path: normalise, demosaic with `method`, apply `OpcodeList3`, then the default crop.
+    /// The result is camera RGB (not white balanced), white level = 1.0, not oriented.
+    pub fn develop(&self, method: Method) -> Result<Rgb32f> {
+        let n = self.normalized()?;
+        let mut rgb = demosaic(&n, method);
+        opcodes::apply_list3(&self.opcodes.list3, &mut rgb);
+        let c = self.crop.clipped(rgb.width, rgb.height);
+        if c.width == 0 || c.height == 0 || (c.x == 0 && c.y == 0 && c.width == rgb.width && c.height == rgb.height) {
+            return Ok(rgb);
+        }
+        Ok(rgb.crop(c.x, c.y, c.width, c.height))
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn cfa_helpers() {
+        let c = Cfa::bayer("RGGB").unwrap();
+        assert!(c.is_bayer());
+        assert_eq!(c.color_at(0, 0), 0);
+        assert_eq!(c.color_at(3, 3), 2);
+        assert_eq!(c.shifted(1, 0).name(), "GRBG");
+        assert_eq!(c.shifted(1, 1).name(), "BGGR");
+        assert!(Cfa::bayer("RGGX").is_none());
+        let x = Cfa::xtrans();
+        assert!(!x.is_bayer());
+        let greens = x.pattern.iter().filter(|&&p| p == 1).count();
+        assert_eq!(greens, 20);
+        // every 3×3 block of the X-Trans layout has 5 greens
+        for by in 0..2 {
+            for bx in 0..2 {
+                let g = (0..9).filter(|i| x.color_at(bx * 3 + i % 3, by * 3 + i / 3) == 1).count();
+                assert_eq!(g, 5);
+            }
+        }
+    }
+
+    #[test]
+    fn black_levels() {
+        let b = BlackLevel { repeat_rows: 2, repeat_cols: 2, values: vec![1.0, 2.0, 3.0, 4.0], delta_h: vec![0.5, 0.0], delta_v: vec![0.0, 10.0] };
+        assert_eq!(b.at(0, 0, 0, 1), 1.5);
+        assert_eq!(b.at(1, 1, 0, 1), 14.0);
+        assert_eq!(b.at(3, 2, 0, 1), 2.0);
+        assert_eq!(b.mean(), 2.5);
+        assert_eq!(BlackLevel::uniform(7.0).at(9, 9, 0, 3), 7.0);
+    }
+
+    #[test]
+    fn probe_rejects_junk() {
+        assert_eq!(probe(b""), None);
+        assert_eq!(probe(b"hello world, not a raw"), None);
+        assert_eq!(probe(b"\0\0\0\x18ftypcrx \0\0\0\x01"), Some(RawFormat::Cr3));
+        assert_eq!(probe(b"FUJIFILMCCD-RAW 0201"), Some(RawFormat::Raf));
+        assert!(matches!(decode(b"FUJIFILMCCD-RAW 0201"), Err(RawError::Unsupported(_))));
+        assert_eq!(decode(b"junk"), Err(RawError::NotRaw));
+    }
+}
diff --git a/crates/raw/src/ljpeg.rs b/crates/raw/src/ljpeg.rs
new file mode 100644
index 00000000..000a2364
--- /dev/null
+++ b/crates/raw/src/ljpeg.rs
@@ -0,0 +1,786 @@
+//! 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`),
+//!   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.
+
+use crate::RawError;
+
+/// A decoded lossless-JPEG frame.
+#[derive(Clone, Debug, PartialEq, Eq)]
+pub struct Frame {
+    pub width: usize,
+    pub height: usize,
+    pub components: usize,
+    pub precision: u8,
+    pub predictor: u8,
+    pub point_transform: u8,
+    /// `width × height × components` interleaved samples.
+    pub data: Vec<u16>,
+}
+
+fn err(s: &str) -> RawError {
+    RawError::Corrupt(format!("lossless JPEG: {s}"))
+}
+
+/// Huffman decoding table (T.81 Annex C "BITS"/"HUFFVAL").
+#[derive(Clone, Debug)]
+pub struct Huffman {
+    /// For peeks of `LOOKUP` bits: (code length, symbol); length 0 = slow path.
+    fast: Vec<(u8, u8)>,
+    maxcode: [i32; 18],
+    valptr: [i32; 17],
+    mincode: [i32; 17],
+    values: Vec<u8>,
+}
+
+const LOOKUP: u32 = 10;
+
+impl Huffman {
+    /// Build from the 16 code-length counts and the symbol list.
+    pub fn new(counts: &[u8; 16], values: &[u8]) -> Result<Huffman, RawError> {
+        let total: usize = counts.iter().map(|&c| c as usize).sum();
+        if total > values.len() || total > 256 || total == 0 {
+            return Err(err("bad Huffman table"));
+        }
+        let mut maxcode = [-1i32; 18];
+        let mut valptr = [0i32; 17];
+        let mut mincode = [0i32; 17];
+        let mut fast = vec![(0u8, 0u8); 1 << LOOKUP];
+        let mut code = 0i32;
+        let mut k = 0usize;
+        for len in 1..=16usize {
+            let n = counts[len - 1] as i32;
+            valptr[len] = k as i32;
+            mincode[len] = code;
+            for _ in 0..n {
+                if len as u32 <= LOOKUP {
+                    let shift = LOOKUP - len as u32;
+                    let base = (code as usize) << shift;
+                    for i in 0..(1usize << shift) {
+                        if let Some(slot) = fast.get_mut(base + i) {
+                            *slot = (len as u8, values[k]);
+                        }
+                    }
+                }
+                code += 1;
+                k += 1;
+            }
+            maxcode[len] = if n > 0 { code - 1 } else { -1 };
+            if code > (1 << len) {
+                return Err(err("over-subscribed Huffman table"));
+            }
+            code <<= 1;
+        }
+        maxcode[17] = i32::MAX;
+        Ok(Huffman { fast, maxcode, valptr, mincode, values: values[..total].to_vec() })
+    }
+
+    #[inline]
+    pub fn decode(&self, br: &mut BitReader) -> Result<u8, RawError> {
+        let peek = br.peek(LOOKUP);
+        let (len, sym) = self.fast[peek as usize];
+        if len > 0 {
+            br.consume(len as u32);
+            return Ok(sym);
+        }
+        let mut code = br.peek(LOOKUP + 1) as i32;
+        let mut len = LOOKUP as usize + 1;
+        while len <= 16 && code > self.maxcode[len] {
+            len += 1;
+            code = br.peek(len as u32) as i32;
+        }
+        if len > 16 {
+            return Err(err("invalid Huffman code"));
+        }
+        br.consume(len as u32);
+        let idx = self.valptr[len] + code - self.mincode[len];
+        self.values.get(idx as usize).copied().ok_or_else(|| err("invalid Huffman code"))
+    }
+}
+
+/// MSB-first bit reader over JPEG entropy-coded data (handles `FF 00` stuffing; stops at markers, feeding zeros).
+pub struct BitReader<'a> {
+    data: &'a [u8],
+    pos: usize,
+    buf: u64,
+    bits: u32,
+    /// A marker (`FF xx`, xx ≠ 0) was reached.
+    pub marker: Option<u8>,
+    /// Bits fed past the end of data (for corruption detection).
+    pub overrun: u32,
+}
+
+impl<'a> BitReader<'a> {
+    pub fn new(data: &'a [u8]) -> Self {
+        Self { data, pos: 0, buf: 0, bits: 0, marker: None, overrun: 0 }
+    }
+    #[inline]
+    fn fill(&mut self) {
+        while self.bits <= 56 {
+            let byte = if self.marker.is_some() || self.pos >= self.data.len() {
+                self.overrun += 8;
+                0
+            } else {
+                let b = self.data[self.pos];
+                if b == 0xff {
+                    let next = self.data.get(self.pos + 1).copied().unwrap_or(0);
+                    if next == 0 {
+                        self.pos += 2;
+                        0xff
+                    } else {
+                        self.marker = Some(next);
+                        self.overrun += 8;
+                        0
+                    }
+                } else {
+                    self.pos += 1;
+                    b
+                }
+            };
+            self.buf |= (byte as u64) << (56 - self.bits);
+            self.bits += 8;
+        }
+    }
+    #[inline]
+    pub fn peek(&mut self, n: u32) -> u32 {
+        if self.bits < n {
+            self.fill();
+        }
+        (self.buf >> (64 - n)) as u32
+    }
+    #[inline]
+    pub fn consume(&mut self, n: u32) {
+        self.buf <<= n;
+        self.bits -= n;
+    }
+    #[inline]
+    pub fn get(&mut self, n: u32) -> u32 {
+        if n == 0 {
+            return 0;
+        }
+        let v = self.peek(n);
+        self.consume(n);
+        v
+    }
+    /// Discard buffered bits and skip to just after the next `RSTn` marker. Returns false if none found.
+    fn restart(&mut self) -> bool {
+        self.buf = 0;
+        self.bits = 0;
+        self.marker = None;
+        self.overrun = 0;
+        while self.pos + 1 < self.data.len() {
+            if self.data[self.pos] == 0xff && (0xd0..=0xd7).contains(&self.data[self.pos + 1]) {
+                self.pos += 2;
+                return true;
+            }
+            self.pos += 1;
+        }
+        false
+    }
+    /// Bytes consumed so far (approximate while bits are buffered).
+    pub fn position(&self) -> usize {
+        self.pos
+    }
+}
+
+/// Difference value for category `ssss` (T.81 Table H.2 / F.12 EXTEND).
+#[inline]
+pub fn diff_value(br: &mut BitReader, ssss: u8) -> i32 {
+    match ssss {
+        0 => 0,
+        16 => 32768,
+        s => {
+            let s = s.min(16) as u32;
+            let v = br.get(s) as i32;
+            if v < (1 << (s - 1)) { v - (1 << s) + 1 } else { v }
+        }
+    }
+}
+
+struct Header {
+    precision: u8,
+    height: usize,
+    width: usize,
+    comps: Vec<u8>,
+    tables: [Option<Huffman>; 4],
+    table_for: Vec<usize>,
+    predictor: u8,
+    pt: u8,
+    restart: usize,
+    scan_start: usize,
+}
+
+fn be16(d: &[u8], i: usize) -> Result<usize, RawError> {
+    d.get(i..i + 2).map(|s| u16::from_be_bytes([s[0], s[1]]) as usize).ok_or_else(|| err("truncated header"))
+}
+
+fn parse_header(d: &[u8]) -> Result<Header, RawError> {
+    if d.len() < 4 || d[0] != 0xff || d[1] != 0xd8 {
+        return Err(err("missing SOI"));
+    }
+    let mut h = Header {
+        precision: 0,
+        height: 0,
+        width: 0,
+        comps: vec![],
+        tables: [None, None, None, None],
+        table_for: vec![],
+        predictor: 1,
+        pt: 0,
+        restart: 0,
+        scan_start: 0,
+    };
+    let mut i = 2;
+    let mut have_sof = false;
+    loop {
+        // skip fill bytes / stray data up to the next marker
+        while i < d.len() && d[i] != 0xff {
+            i += 1;
+        }
+        while i < d.len() && d[i] == 0xff {
+            i += 1;
+        }
+        let Some(&m) = d.get(i) else { return Err(err("no scan")) };
+        i += 1;
+        if m == 0xd8 || (0xd0..=0xd7).contains(&m) {
+            continue;
+        }
+        if m == 0xd9 {
+            return Err(err("EOI before scan"));
+        }
+        let len = be16(d, i)?;
+        if len < 2 || i + len > d.len() {
+            return Err(err("truncated segment"));
+        }
+        let seg = &d[i + 2..i + len];
+        match m {
+            0xc3 => {
+                if seg.len() < 6 {
+                    return Err(err("short SOF3"));
+                }
+                h.precision = seg[0];
+                h.height = u16::from_be_bytes([seg[1], seg[2]]) as usize;
+                h.width = u16::from_be_bytes([seg[3], seg[4]]) as usize;
+                let n = seg[5] as usize;
+                if n == 0 || n > 4 || seg.len() < 6 + 3 * n {
+                    return Err(err("bad component count"));
+                }
+                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()));
+                    }
+                    h.comps.push(id);
+                }
+                if !(2..=16).contains(&h.precision) {
+                    return Err(err("bad precision"));
+                }
+                have_sof = true;
+            }
+            0xc0..=0xc2 | 0xc5..=0xc7 | 0xc9..=0xcb | 0xcd..=0xcf => {
+                return Err(RawError::Unsupported(format!("JPEG process SOF{:X} (not lossless Huffman)", m - 0xc0)));
+            }
+            0xc4 => {
+                let mut j = 0;
+                while j + 17 <= seg.len() {
+                    let (class, id) = (seg[j] >> 4, (seg[j] & 15) as usize);
+                    let mut counts = [0u8; 16];
+                    counts.copy_from_slice(&seg[j + 1..j + 17]);
+                    let n: usize = counts.iter().map(|&c| c as usize).sum();
+                    let vals = seg.get(j + 17..j + 17 + n).ok_or_else(|| err("truncated DHT"))?;
+                    if class == 0 && id < 4 {
+                        h.tables[id] = Some(Huffman::new(&counts, vals)?);
+                    }
+                    j += 17 + n;
+                }
+            }
+            0xdd => h.restart = be16(seg, 0)?,
+            0xda => {
+                if !have_sof {
+                    return Err(err("SOS before SOF3"));
+                }
+                let ns = *seg.first().ok_or_else(|| err("short SOS"))? as usize;
+                if ns != h.comps.len() || seg.len() < 1 + 2 * ns + 3 {
+                    return Err(RawError::Unsupported("lossless JPEG with multiple scans".into()));
+                }
+                for c in 0..ns {
+                    let cid = seg[1 + 2 * c];
+                    let pos = h.comps.iter().position(|&x| x == cid).ok_or_else(|| err("unknown scan component"))?;
+                    if pos != c {
+                        return Err(err("scan component order"));
+                    }
+                    let td = (seg[2 + 2 * c] >> 4) as usize;
+                    if td > 3 || h.tables[td].is_none() {
+                        return Err(err("missing Huffman table"));
+                    }
+                    h.table_for.push(td);
+                }
+                h.predictor = seg[1 + 2 * ns];
+                h.pt = seg[3 + 2 * ns] & 15;
+                if !(0..=7).contains(&h.predictor) {
+                    return Err(err("bad predictor"));
+                }
+                if h.pt >= h.precision {
+                    return Err(err("bad point transform"));
+                }
+                h.scan_start = i + len;
+                return Ok(h);
+            }
+            _ => {}
+        }
+        i += len;
+    }
+}
+
+/// Parse only the frame header: (width, height, components, precision).
+pub fn frame_info(d: &[u8]) -> Result<(usize, usize, usize, u8), RawError> {
+    let h = parse_header(d)?;
+    Ok((h.width, h.height, h.comps.len(), h.precision))
+}
+
+/// Decode a lossless JPEG stream. `max_samples` bounds the allocation (hostile headers).
+pub fn decode(d: &[u8], max_samples: usize) -> Result<Frame, 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 total = w.checked_mul(ht).and_then(|v| v.checked_mul(nc)).ok_or_else(|| err("frame too large"))?;
+    if total > max_samples {
+        return Err(RawError::Limit("lossless JPEG frame larger than expected"));
+    }
+    // The entropy data must hold at least ~1 bit per sample; reject absurd headers early.
+    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"));
+    }
+    let tables: Vec<&Huffman> = h.table_for.iter().map(|&t| h.tables[t].as_ref().expect("checked")).collect();
+    let mut out = vec![0u16; total];
+    let mut br = BitReader::new(entropy);
+    let init = 1i32 << (h.precision - h.pt - 1);
+    let mask = 0xffffi32;
+    let row_len = w * nc;
+    let pred = h.predictor;
+    let mut restart_left = h.restart;
+    let mut first_line = true; // first line of image or restart interval
+    let mut reset_pending = true;
+    for y in 0..ht {
+        let (prev_rows, cur_rows) = out.split_at_mut(y * row_len);
+        let cur = &mut cur_rows[..row_len];
+        let prev = if y > 0 { &prev_rows[(y - 1) * row_len..] } else { &[][..] };
+        for x in 0..w {
+            if h.restart > 0 {
+                if restart_left == 0 {
+                    if !br.restart() {
+                        return Err(err("missing restart marker"));
+                    }
+                    restart_left = h.restart;
+                    first_line = true;
+                    reset_pending = true;
+                }
+                restart_left -= 1;
+            }
+            for c in 0..nc {
+                let i = x * nc + c;
+                let p = if reset_pending {
+                    init
+                } else if first_line {
+                    if x == 0 { prev.get(i).map_or(init, |&v| v as i32) } else { cur[i - nc] as i32 }
+                } else if x == 0 || pred == 0 {
+                    if pred == 0 { 0 } else { prev[i] as i32 }
+                } else {
+                    let ra = cur[i - nc] as i32;
+                    let rb = prev[i] as i32;
+                    let rc = prev[i - nc] 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);
+                cur[i] = ((p + diff) & mask) as u16;
+            }
+            reset_pending = false;
+        }
+        first_line = false;
+        if br.overrun > 64 * 8 {
+            return Err(err("entropy data exhausted"));
+        }
+    }
+    if h.pt > 0 {
+        for v in &mut out {
+            *v <<= h.pt;
+        }
+    }
+    Ok(Frame { width: w, height: ht, components: nc, precision: h.precision, predictor: pred, point_transform: h.pt, data: out })
+}
+
+// ------------------------------------------------------------------------------------------------ encoder
+
+/// Build Huffman code lengths from frequencies (T.81 Annex K.2, Figures K.1–K.3), limited to 16 bits.
+fn code_lengths(freq_in: &[u32; 17]) -> ([u8; 16], Vec<u8>) {
+    let mut freq = [0i64; 18];
+    for (i, &f) in freq_in.iter().enumerate() {
+        freq[i] = f as i64;
+    }
+    freq[17] = 1; // reserved symbol so no code is all ones
+    let mut codesize = [0usize; 18];
+    let mut others = [-1i32; 18];
+    loop {
+        // V1: least frequency > 0 (largest index on ties), V2: next least
+        let mut v1 = -1i32;
+        let mut v2 = -1i32;
+        for i in 0..18 {
+            if freq[i] > 0 && (v1 < 0 || freq[i] <= freq[v1 as usize]) {
+                v1 = i as i32;
+            }
+        }
+        for i in 0..18 {
+            if freq[i] > 0 && i as i32 != v1 && (v2 < 0 || freq[i] <= freq[v2 as usize]) {
+                v2 = i as i32;
+            }
+        }
+        if v2 < 0 {
+            break;
+        }
+        let (a, b) = (v1 as usize, v2 as usize);
+        freq[a] += freq[b];
+        freq[b] = 0;
+        let mut x = a;
+        codesize[x] += 1;
+        while others[x] >= 0 {
+            x = others[x] as usize;
+            codesize[x] += 1;
+        }
+        others[x] = b as i32;
+        let mut x = b;
+        codesize[x] += 1;
+        while others[x] >= 0 {
+            x = others[x] as usize;
+            codesize[x] += 1;
+        }
+    }
+    let mut bits = [0i32; 33];
+    for &cs in &codesize {
+        if cs > 0 {
+            bits[cs.min(32)] += 1;
+        }
+    }
+    // Adjust_BITS (Figure K.3)
+    let mut i = 32;
+    while i > 16 {
+        while bits[i] > 0 {
+            let mut j = i - 2;
+            while bits[j] == 0 {
+                j -= 1;
+            }
+            bits[i] -= 2;
+            bits[i - 1] += 1;
+            bits[j + 1] += 2;
+            bits[j] -= 1;
+        }
+        i -= 1;
+    }
+    while bits[i] == 0 {
+        i -= 1;
+    }
+    bits[i] -= 1; // remove the reserved code
+    let mut counts = [0u8; 16];
+    for l in 1..=16 {
+        counts[l - 1] = bits[l] as u8;
+    }
+    // Sort_input (Figure K.4): symbols by code size
+    let mut values = Vec::new();
+    for l in 1..=32 {
+        for (s, &cs) in codesize.iter().enumerate().take(17) {
+            if cs == l {
+                values.push(s as u8);
+            }
+        }
+    }
+    (counts, values)
+}
+
+/// Canonical codes for a table: symbol → (code, length).
+fn canonical_codes(counts: &[u8; 16], values: &[u8]) -> [(u32, u8); 17] {
+    let mut table = [(0u32, 0u8); 17];
+    let mut code = 0u32;
+    let mut k = 0;
+    for len in 1..=16u8 {
+        for _ in 0..counts[len as usize - 1] {
+            if let Some(&s) = values.get(k)
+                && (s as usize) < 17
+            {
+                table[s as usize] = (code, len);
+            }
+            code += 1;
+            k += 1;
+        }
+        code <<= 1;
+    }
+    table
+}
+
+struct BitWriter {
+    out: Vec<u8>,
+    acc: u64,
+    n: u32,
+}
+
+impl BitWriter {
+    fn put(&mut self, v: u32, bits: u32) {
+        if bits == 0 {
+            return;
+        }
+        self.acc = (self.acc << bits) | (v as u64 & ((1u64 << bits) - 1));
+        self.n += bits;
+        while self.n >= 8 {
+            let b = (self.acc >> (self.n - 8)) as u8;
+            self.out.push(b);
+            if b == 0xff {
+                self.out.push(0);
+            }
+            self.n -= 8;
+        }
+    }
+    fn flush(&mut self) {
+        if self.n > 0 {
+            let pad = 8 - self.n;
+            self.put((1 << pad) - 1, pad);
+        }
+    }
+}
+
+fn ssss_of(diff: i32) -> u8 {
+    if diff == 0 {
+        0
+    } else if diff == 32768 || diff == -32768 {
+        16
+    } else {
+        32 - (diff.unsigned_abs()).leading_zeros() as u8
+    }
+}
+
+/// Encode samples (`width × height × components`, interleaved) as a lossless JPEG with `predictor` (1–7) and
+/// `precision` bits. `restart` > 0 inserts restart markers every `restart` MCUs (pixels).
+pub fn encode(data: &[u16], width: usize, height: usize, components: usize, precision: u8, predictor: u8, restart: usize) -> Vec<u8> {
+    assert!((1..=4).contains(&components) && (2..=16).contains(&precision) && (1..=7).contains(&predictor));
+    assert_eq!(data.len(), width * height * components);
+    let init = 1i32 << (precision - 1);
+    let row_len = width * components;
+    // compute differences in the same order as the decoder
+    let mut diffs = vec![0i32; data.len()];
+    let mut restart_left = restart;
+    let mut first_line = true;
+    let mut reset_pending = true;
+    let mut marks = Vec::new(); // sample index where an RST marker precedes
+    for y in 0..height {
+        for x in 0..width {
+            if restart > 0 {
+                if restart_left == 0 {
+                    marks.push((y * width + x) * components);
+                    restart_left = restart;
+                    first_line = true;
+                    reset_pending = true;
+                }
+                restart_left -= 1;
+            }
+            for c in 0..components {
+                let i = y * row_len + x * components + c;
+                let p = if reset_pending {
+                    init
+                } else if first_line {
+                    if x == 0 { data[i - row_len] as i32 } else { data[i - components] as i32 }
+                } else if x == 0 {
+                    data[i - row_len] as i32
+                } else {
+                    let ra = data[i - components] as i32;
+                    let rb = data[i - row_len] as i32;
+                    let rc = data[i - row_len - components] 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,
+                    }
+                };
+                // modulo 2^16 difference mapped into [-32767, 32768]
+                let mut d = (data[i] as i32 - p) & 0xffff;
+                if d > 32768 {
+                    d -= 65536;
+                }
+                diffs[i] = d;
+            }
+            reset_pending = false;
+        }
+        first_line = false;
+    }
+    let mut tables = Vec::new();
+    for c in 0..components {
+        let mut freq = [0u32; 17];
+        for d in diffs.iter().skip(c).step_by(components) {
+            freq[ssss_of(*d) as usize] += 1;
+        }
+        let (counts, values) = code_lengths(&freq);
+        tables.push((counts, values));
+    }
+    let mut out = vec![0xff, 0xd8];
+    // SOF3
+    out.extend_from_slice(&[0xff, 0xc3]);
+    out.extend_from_slice(&((8 + 3 * components) as u16).to_be_bytes());
+    out.push(precision);
+    out.extend_from_slice(&(height as u16).to_be_bytes());
+    out.extend_from_slice(&(width as u16).to_be_bytes());
+    out.push(components as u8);
+    for c in 0..components {
+        out.extend_from_slice(&[c as u8 + 1, 0x11, 0]);
+    }
+    // DHT
+    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());
+    }
+    // SOS
+    out.extend_from_slice(&[0xff, 0xda]);
+    out.extend_from_slice(&((6 + 2 * components) as u16).to_be_bytes());
+    out.push(components as u8);
+    for c in 0..components {
+        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 mark_iter = marks.iter().peekable();
+    let mut rst = 0u8;
+    for (i, &d) in diffs.iter().enumerate() {
+        if mark_iter.peek().is_some_and(|&&m| m == i) {
+            mark_iter.next();
+            bw.flush();
+            bw.out.extend_from_slice(&[0xff, 0xd0 + rst]);
+            rst = (rst + 1) & 7;
+            bw.acc = 0;
+        }
+        let c = i % components;
+        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::*;
+
+    fn noise(n: usize, bits: u32, seed: u64) -> Vec<u16> {
+        let mut s = seed;
+        (0..n)
+            .map(|i| {
+                s = s.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407);
+                let smooth = ((i % 97) as u32 * 37) % (1 << bits);
+                let jitter = ((s >> 33) as u32) % 64;
+                ((smooth + jitter) % (1 << bits)) as u16
+            })
+            .collect()
+    }
+
+    #[test]
+    fn roundtrip_all_predictors_components_precisions() {
+        for pred in 1..=7u8 {
+            for comps in 1..=4usize {
+                for bits in [8u8, 12, 14, 16] {
+                    let (w, h) = (13, 7);
+                    let data = noise(w * h * comps, bits as u32, (pred as u64) * 31 + comps as u64);
+                    let enc = encode(&data, w, h, comps, bits, pred, 0);
+                    let f = decode(&enc, usize::MAX).unwrap();
+                    assert_eq!((f.width, f.height, f.components, f.precision, f.predictor), (w, h, comps, bits, pred));
+                    assert_eq!(f.data, data, "pred {pred} comps {comps} bits {bits}");
+                }
+            }
+        }
+    }
+
+    #[test]
+    fn roundtrip_extremes_and_restarts() {
+        // alternating 0 / 65535 forces SSSS = 16 differences
+        let data: Vec<u16> = (0..64).map(|i| if i % 2 == 0 { 0 } else { 65535 }).collect();
+        for pred in 1..=7 {
+            let enc = encode(&data, 8, 8, 1, 16, pred, 0);
+            assert_eq!(decode(&enc, usize::MAX).unwrap().data, data);
+        }
+        let data = noise(20 * 9 * 2, 12, 5);
+        for rst in [1usize, 7, 20, 40] {
+            let enc = encode(&data, 20, 9, 2, 12, 6, rst);
+            assert_eq!(decode(&enc, usize::MAX).unwrap().data, data, "restart {rst}");
+        }
+    }
+
+    #[test]
+    fn constant_image_single_symbol_table() {
+        let data = vec![1234u16; 16 * 4];
+        let enc = encode(&data, 16, 4, 1, 14, 1, 0);
+        assert_eq!(decode(&enc, usize::MAX).unwrap().data, data);
+    }
+
+    #[test]
+    fn rejects_bad_streams() {
+        assert!(decode(&[], 100).is_err());
+        assert!(decode(&[0xff, 0xd8, 0xff, 0xd9], 100).is_err());
+        let data = noise(64, 12, 1);
+        let enc = encode(&data, 8, 8, 1, 12, 1, 0);
+        assert!(matches!(decode(&enc, 10), Err(RawError::Limit(_))));
+        for n in 0..enc.len() {
+            let _ = decode(&enc[..n], usize::MAX);
+        }
+        let mut bad = enc.clone();
+        for i in (20..bad.len()).step_by(3) {
+            bad[i] ^= 0x5a;
+        }
+        let _ = decode(&bad, usize::MAX);
+    }
+
+    #[test]
+    fn code_lengths_are_limited_to_16() {
+        // Fibonacci-like frequencies produce very long codes before adjustment
+        let mut f = [0u32; 17];
+        let (mut a, mut b) = (1u32, 1u32);
+        for x in f.iter_mut() {
+            *x = a;
+            let c = a.saturating_add(b);
+            a = b;
+            b = c;
+        }
+        let (counts, values) = code_lengths(&f);
+        assert_eq!(counts.iter().map(|&c| c as usize).sum::<usize>(), 17);
+        assert_eq!(values.len(), 17);
+        assert!(Huffman::new(&counts, &values).is_ok());
+    }
+}
diff --git a/crates/raw/src/opcodes.rs b/crates/raw/src/opcodes.rs
new file mode 100644
index 00000000..03ceebc8
--- /dev/null
+++ b/crates/raw/src/opcodes.rs
@@ -0,0 +1,713 @@
+//! DNG opcode lists (`OpcodeList1/2/3`, DNG 1.7 chapter 7). Opcode lists are always big-endian.
+//!
+//! Parsed: every opcode (unknown ones are kept as [`Opcode::Unknown`]). Applied: `WarpRectilinear`,
+//! `FixVignetteRadial`, `FixBadPixelsConstant`, `FixBadPixelsList`, `MapTable`, `MapPolynomial`, `GainMap`,
+//! `DeltaPerRow/Column`, `ScalePerRow/Column`. `TrimBounds` and `WarpFisheye` are recorded but not applied.
+//!
+//! Value convention: list 1 runs on raw sample values (16-bit scale), lists 2 and 3 on values normalised to
+//! [0, 1]; table/polynomial/delta opcodes are defined on the normalised range and are rescaled accordingly.
+
+use crate::{Cfa, Rgb32f};
+use serde::{Deserialize, Serialize};
+
+/// Rectangular area + plane selection + pitch shared by the "area" opcodes.
+#[derive(Clone, Copy, Debug, Default, PartialEq, Eq, Serialize, Deserialize)]
+pub struct Area {
+    pub top: u32,
+    pub left: u32,
+    pub bottom: u32,
+    pub right: u32,
+    pub plane: u32,
+    pub planes: u32,
+    pub row_pitch: u32,
+    pub col_pitch: u32,
+}
+
+#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
+pub enum Opcode {
+    /// Per-plane radial (kr0..kr3) + tangential (kt0, kt1) coefficients and the relative optical centre.
+    WarpRectilinear {
+        planes: Vec<[f64; 6]>,
+        center: [f64; 2],
+    },
+    WarpFisheye {
+        planes: Vec<[f64; 4]>,
+        center: [f64; 2],
+    },
+    FixVignetteRadial {
+        k: [f64; 5],
+        center: [f64; 2],
+    },
+    FixBadPixelsConstant {
+        constant: u32,
+        bayer_phase: u32,
+    },
+    FixBadPixelsList {
+        bayer_phase: u32,
+        points: Vec<(u32, u32)>,
+        rects: Vec<[u32; 4]>,
+    },
+    TrimBounds {
+        top: u32,
+        left: u32,
+        bottom: u32,
+        right: u32,
+    },
+    MapTable {
+        area: Area,
+        table: Vec<u16>,
+    },
+    MapPolynomial {
+        area: Area,
+        coefficients: Vec<f64>,
+    },
+    GainMap {
+        area: Area,
+        points_v: u32,
+        points_h: u32,
+        spacing: [f64; 2],
+        origin: [f64; 2],
+        map_planes: u32,
+        gains: Vec<f32>,
+    },
+    DeltaPerRow {
+        area: Area,
+        deltas: Vec<f32>,
+    },
+    DeltaPerColumn {
+        area: Area,
+        deltas: Vec<f32>,
+    },
+    ScalePerRow {
+        area: Area,
+        scales: Vec<f32>,
+    },
+    ScalePerColumn {
+        area: Area,
+        scales: Vec<f32>,
+    },
+    Unknown {
+        id: u32,
+        flags: u32,
+        params: Vec<u8>,
+    },
+}
+
+impl Opcode {
+    /// Whether [`apply_list`] / [`apply_list3`] implement this opcode.
+    pub fn is_applied(&self) -> bool {
+        !matches!(self, Opcode::Unknown { .. } | Opcode::TrimBounds { .. } | Opcode::WarpFisheye { .. })
+    }
+}
+
+/// The three opcode lists of a DNG raw IFD.
+#[derive(Clone, Debug, Default, PartialEq, Serialize, Deserialize)]
+pub struct OpcodeLists {
+    pub list1: Vec<Opcode>,
+    pub list2: Vec<Opcode>,
+    pub list3: Vec<Opcode>,
+}
+
+struct Rd<'a> {
+    d: &'a [u8],
+    p: usize,
+}
+
+impl Rd<'_> {
+    fn u32(&mut self) -> Option<u32> {
+        let v = u32::from_be_bytes(self.d.get(self.p..self.p + 4)?.try_into().ok()?);
+        self.p += 4;
+        Some(v)
+    }
+    fn u16(&mut self) -> Option<u16> {
+        let v = u16::from_be_bytes(self.d.get(self.p..self.p + 2)?.try_into().ok()?);
+        self.p += 2;
+        Some(v)
+    }
+    fn f64(&mut self) -> Option<f64> {
+        let v = f64::from_be_bytes(self.d.get(self.p..self.p + 8)?.try_into().ok()?);
+        self.p += 8;
+        Some(v)
+    }
+    fn f32(&mut self) -> Option<f32> {
+        let v = f32::from_be_bytes(self.d.get(self.p..self.p + 4)?.try_into().ok()?);
+        self.p += 4;
+        Some(v)
+    }
+    fn area(&mut self) -> Option<Area> {
+        Some(Area {
+            top: self.u32()?,
+            left: self.u32()?,
+            bottom: self.u32()?,
+            right: self.u32()?,
+            plane: self.u32()?,
+            planes: self.u32()?,
+            row_pitch: self.u32()?.max(1),
+            col_pitch: self.u32()?.max(1),
+        })
+    }
+    fn remaining(&self) -> usize {
+        self.d.len().saturating_sub(self.p)
+    }
+}
+
+fn parse_one(id: u32, flags: u32, params: &[u8]) -> Option<Opcode> {
+    let mut r = Rd { d: params, p: 0 };
+    let f = |r: &mut Rd| r.f64().filter(|v| v.is_finite());
+    Some(match id {
+        1 => {
+            let n = r.u32()? as usize;
+            if n == 0 || n > 4 {
+                return None;
+            }
+            let mut planes = Vec::with_capacity(n);
+            for _ in 0..n {
+                planes.push([f(&mut r)?, f(&mut r)?, f(&mut r)?, f(&mut r)?, f(&mut r)?, f(&mut r)?]);
+            }
+            Opcode::WarpRectilinear { planes, center: [f(&mut r)?, f(&mut r)?] }
+        }
+        2 => {
+            let n = r.u32()? as usize;
+            if n == 0 || n > 4 {
+                return None;
+            }
+            let mut planes = Vec::with_capacity(n);
+            for _ in 0..n {
+                planes.push([f(&mut r)?, f(&mut r)?, f(&mut r)?, f(&mut r)?]);
+            }
+            Opcode::WarpFisheye { planes, center: [f(&mut r)?, f(&mut r)?] }
+        }
+        3 => Opcode::FixVignetteRadial { k: [f(&mut r)?, f(&mut r)?, f(&mut r)?, f(&mut r)?, f(&mut r)?], center: [f(&mut r)?, f(&mut r)?] },
+        4 => Opcode::FixBadPixelsConstant { constant: r.u32()?, bayer_phase: r.u32()? },
+        5 => {
+            let bayer_phase = r.u32()?;
+            let (np, nr) = (r.u32()? as usize, r.u32()? as usize);
+            if np.saturating_mul(8).saturating_add(nr.saturating_mul(16)) > r.remaining() {
+                return None;
+            }
+            let points = (0..np).map(|_| Some((r.u32()?, r.u32()?))).collect::<Option<_>>()?;
+            let rects = (0..nr).map(|_| Some([r.u32()?, r.u32()?, r.u32()?, r.u32()?])).collect::<Option<_>>()?;
+            Opcode::FixBadPixelsList { bayer_phase, points, rects }
+        }
+        6 => Opcode::TrimBounds { top: r.u32()?, left: r.u32()?, bottom: r.u32()?, right: r.u32()? },
+        7 => {
+            let area = r.area()?;
+            let n = r.u32()? as usize;
+            if n == 0 || n > 65536 || n * 2 > r.remaining() {
+                return None;
+            }
+            Opcode::MapTable { area, table: (0..n).map(|_| r.u16()).collect::<Option<_>>()? }
+        }
+        8 => {
+            let area = r.area()?;
+            let deg = r.u32()? as usize;
+            if deg > 8 {
+                return None;
+            }
+            Opcode::MapPolynomial { area, coefficients: (0..=deg).map(|_| f(&mut r)).collect::<Option<_>>()? }
+        }
+        9 => {
+            let area = r.area()?;
+            let (pv, ph) = (r.u32()?, r.u32()?);
+            let spacing = [f(&mut r)?, f(&mut r)?];
+            let origin = [f(&mut r)?, f(&mut r)?];
+            let mp = r.u32()?;
+            let n = (pv as usize).checked_mul(ph as usize)?.checked_mul(mp as usize)?;
+            if n == 0 || n * 4 > r.remaining() {
+                return None;
+            }
+            let gains = (0..n).map(|_| r.f32()).collect::<Option<_>>()?;
+            Opcode::GainMap { area, points_v: pv, points_h: ph, spacing, origin, map_planes: mp, gains }
+        }
+        10..=13 => {
+            let area = r.area()?;
+            let n = r.u32()? as usize;
+            if n * 4 > r.remaining() {
+                return None;
+            }
+            let v: Vec<f32> = (0..n).map(|_| r.f32()).collect::<Option<_>>()?;
+            match id {
+                10 => Opcode::DeltaPerRow { area, deltas: v },
+                11 => Opcode::DeltaPerColumn { area, deltas: v },
+                12 => Opcode::ScalePerRow { area, scales: v },
+                _ => Opcode::ScalePerColumn { area, scales: v },
+            }
+        }
+        _ => Opcode::Unknown { id, flags, params: params.to_vec() },
+    })
+}
+
+fn put_area(o: &mut Vec<u8>, a: &Area) {
+    for v in [a.top, a.left, a.bottom, a.right, a.plane, a.planes, a.row_pitch, a.col_pitch] {
+        o.extend_from_slice(&v.to_be_bytes());
+    }
+}
+
+fn params(op: &Opcode) -> (u32, u32, Vec<u8>) {
+    let mut o = Vec::new();
+    let f64s = |o: &mut Vec<u8>, v: &[f64]| v.iter().for_each(|x| o.extend_from_slice(&x.to_be_bytes()));
+    let u32s = |o: &mut Vec<u8>, v: &[u32]| v.iter().for_each(|x| o.extend_from_slice(&x.to_be_bytes()));
+    let f32s = |o: &mut Vec<u8>, v: &[f32]| v.iter().for_each(|x| o.extend_from_slice(&x.to_be_bytes()));
+    let id = match op {
+        Opcode::WarpRectilinear { planes, center } => {
+            u32s(&mut o, &[planes.len() as u32]);
+            planes.iter().for_each(|p| f64s(&mut o, p));
+            f64s(&mut o, center);
+            1
+        }
+        Opcode::WarpFisheye { planes, center } => {
+            u32s(&mut o, &[planes.len() as u32]);
+            planes.iter().for_each(|p| f64s(&mut o, p));
+            f64s(&mut o, center);
+            2
+        }
+        Opcode::FixVignetteRadial { k, center } => {
+            f64s(&mut o, k);
+            f64s(&mut o, center);
+            3
+        }
+        Opcode::FixBadPixelsConstant { constant, bayer_phase } => {
+            u32s(&mut o, &[*constant, *bayer_phase]);
+            4
+        }
+        Opcode::FixBadPixelsList { bayer_phase, points, rects } => {
+            u32s(&mut o, &[*bayer_phase, points.len() as u32, rects.len() as u32]);
+            points.iter().for_each(|&(r, c)| u32s(&mut o, &[r, c]));
+            rects.iter().for_each(|r| u32s(&mut o, r));
+            5
+        }
+        Opcode::TrimBounds { top, left, bottom, right } => {
+            u32s(&mut o, &[*top, *left, *bottom, *right]);
+            6
+        }
+        Opcode::MapTable { area, table } => {
+            put_area(&mut o, area);
+            u32s(&mut o, &[table.len() as u32]);
+            table.iter().for_each(|v| o.extend_from_slice(&v.to_be_bytes()));
+            7
+        }
+        Opcode::MapPolynomial { area, coefficients } => {
+            put_area(&mut o, area);
+            u32s(&mut o, &[coefficients.len().saturating_sub(1) as u32]);
+            f64s(&mut o, coefficients);
+            8
+        }
+        Opcode::GainMap { area, points_v, points_h, spacing, origin, map_planes, gains } => {
+            put_area(&mut o, area);
+            u32s(&mut o, &[*points_v, *points_h]);
+            f64s(&mut o, spacing);
+            f64s(&mut o, origin);
+            u32s(&mut o, &[*map_planes]);
+            f32s(&mut o, gains);
+            9
+        }
+        Opcode::DeltaPerRow { area, deltas: v }
+        | Opcode::DeltaPerColumn { area, deltas: v }
+        | Opcode::ScalePerRow { area, scales: v }
+        | Opcode::ScalePerColumn { area, scales: v } => {
+            put_area(&mut o, area);
+            u32s(&mut o, &[v.len() as u32]);
+            f32s(&mut o, v);
+            match op {
+                Opcode::DeltaPerRow { .. } => 10,
+                Opcode::DeltaPerColumn { .. } => 11,
+                Opcode::ScalePerRow { .. } => 12,
+                _ => 13,
+            }
+        }
+        Opcode::Unknown { id, flags, params } => return (*id, *flags, params.clone()),
+    };
+    (id, 0, o)
+}
+
+/// Serialise an opcode list (inverse of [`parse_list`]).
+pub fn write_list(list: &[Opcode]) -> Vec<u8> {
+    let mut out = (list.len() as u32).to_be_bytes().to_vec();
+    for op in list {
+        let (id, flags, p) = params(op);
+        for v in [id, 0x0103_0000, flags, p.len() as u32] {
+            out.extend_from_slice(&v.to_be_bytes());
+        }
+        out.extend_from_slice(&p);
+    }
+    out
+}
+
+/// Parse an opcode list blob. Malformed opcodes become [`Opcode::Unknown`]; a truncated list stops early.
+pub fn parse_list(d: &[u8]) -> Vec<Opcode> {
+    let mut r = Rd { d, p: 0 };
+    let Some(count) = r.u32() else { return vec![] };
+    let mut out = Vec::new();
+    for _ in 0..count.min(4096) {
+        let (Some(id), Some(_ver), Some(flags), Some(size)) = (r.u32(), r.u32(), r.u32(), r.u32()) else { break };
+        let Some(params) = d.get(r.p..r.p.saturating_add(size as usize)) else { break };
+        r.p += size as usize;
+        out.push(parse_one(id, flags, params).unwrap_or(Opcode::Unknown { id, flags, params: params.to_vec() }));
+    }
+    out
+}
+
+fn area_iter(a: &Area, w: usize, h: usize, cpp: usize) -> impl Iterator<Item = (usize, usize, usize)> + '_ {
+    let (t, l) = (a.top as usize, a.left as usize);
+    let (b, r) = ((a.bottom as usize).min(h), (a.right as usize).min(w));
+    let (p0, p1) = (a.plane as usize, (a.plane as usize + a.planes.max(1) as usize).min(cpp));
+    (t..b.max(t))
+        .step_by(a.row_pitch.max(1) as usize)
+        .flat_map(move |y| (l..r.max(l)).step_by(a.col_pitch.max(1) as usize).flat_map(move |x| (p0..p1.max(p0)).map(move |p| (x, y, p))))
+}
+
+fn gain_at(points_v: u32, points_h: u32, spacing: [f64; 2], origin: [f64; 2], map_planes: u32, gains: &[f32], rv: f64, rh: f64, plane: usize) -> f32 {
+    let (pv, ph) = (points_v as usize, points_h as usize);
+    let mp = map_planes as usize;
+    let p = plane.min(mp - 1);
+    let fy = if spacing[0] > 0.0 { ((rv - origin[0]) / spacing[0]).clamp(0.0, (pv - 1) as f64) } else { 0.0 };
+    let fx = if spacing[1] > 0.0 { ((rh - origin[1]) / spacing[1]).clamp(0.0, (ph - 1) as f64) } else { 0.0 };
+    let (y0, x0) = (fy.floor() as usize, fx.floor() as usize);
+    let (y1, x1) = ((y0 + 1).min(pv - 1), (x0 + 1).min(ph - 1));
+    let (ty, tx) = ((fy - y0 as f64) as f32, (fx - x0 as f64) as f32);
+    let g = |y: usize, x: usize| gains[(y * ph + x) * mp + p];
+    let top = g(y0, x0) + (g(y0, x1) - g(y0, x0)) * tx;
+    let bot = g(y1, x0) + (g(y1, x1) - g(y1, x0)) * tx;
+    top + (bot - top) * ty
+}
+
+/// Radius normaliser: distance from the centre to the farthest image corner.
+fn max_radius(cx: f64, cy: f64, w: f64, h: f64) -> f64 {
+    [(0.0, 0.0), (w, 0.0), (0.0, h), (w, h)].iter().map(|&(x, y)| (x - cx).hypot(y - cy)).fold(0.0, f64::max).max(1e-9)
+}
+
+/// Apply a list to an interleaved `w × h × cpp` buffer. `cfa` (when the data is a mosaic) is used by the
+/// bad-pixel opcodes. `scale` is the value of 1.0 in the buffer: 65535 for list 1 (raw 16-bit values), 1 for
+/// lists 2 and 3 (normalised values).
+pub fn apply_list(list: &[Opcode], buf: &mut [f32], w: usize, h: usize, cpp: usize, cfa: Option<&Cfa>, scale: f32) {
+    for op in list {
+        match op {
+            Opcode::MapTable { area, table } => {
+                let n = table.len();
+                for (x, y, p) in area_iter(area, w, h, cpp) {
+                    let v = &mut buf[(y * w + x) * cpp + p];
+                    let idx = ((*v / scale * 65535.0).round().max(0.0) as usize).min(n - 1);
+                    *v = table[idx] as f32 / 65535.0 * scale;
+                }
+            }
+            Opcode::MapPolynomial { area, coefficients } => {
+                for (x, y, p) in area_iter(area, w, h, cpp) {
+                    let v = &mut buf[(y * w + x) * cpp + p];
+                    let t = (*v / scale) as f64;
+                    let r = coefficients.iter().rev().fold(0.0, |acc, &c| acc * t + c);
+                    *v = r as f32 * scale;
+                }
+            }
+            Opcode::GainMap { area, points_v, points_h, spacing, origin, map_planes, gains } => {
+                for (x, y, p) in area_iter(area, w, h, cpp) {
+                    let rv = y as f64 / h as f64;
+                    let rh = x as f64 / w as f64;
+                    let g = gain_at(*points_v, *points_h, *spacing, *origin, *map_planes, gains, rv, rh, p - area.plane as usize);
+                    buf[(y * w + x) * cpp + p] *= g;
+                }
+            }
+            Opcode::DeltaPerRow { area, deltas } | Opcode::DeltaPerColumn { area, deltas } => {
+                let per_row = matches!(op, Opcode::DeltaPerRow { .. });
+                for (x, y, p) in area_iter(area, w, h, cpp) {
+                    let i =
+                        if per_row { (y - area.top as usize) / area.row_pitch as usize } else { (x - area.left as usize) / area.col_pitch as usize };
+                    if let Some(d) = deltas.get(i) {
+                        buf[(y * w + x) * cpp + p] += d * scale;
+                    }
+                }
+            }
+            Opcode::ScalePerRow { area, scales } | Opcode::ScalePerColumn { area, scales } => {
+                let per_row = matches!(op, Opcode::ScalePerRow { .. });
+                for (x, y, p) in area_iter(area, w, h, cpp) {
+                    let i =
+                        if per_row { (y - area.top as usize) / area.row_pitch as usize } else { (x - area.left as usize) / area.col_pitch as usize };
+                    if let Some(s) = scales.get(i) {
+                        buf[(y * w + x) * cpp + p] *= s;
+                    }
+                }
+            }
+            Opcode::FixVignetteRadial { k, center } => {
+                let (cx, cy) = (center[0] * w as f64, center[1] * h as f64);
+                let m = max_radius(cx, cy, w as f64, h as f64);
+                for y in 0..h {
+                    for x in 0..w {
+                        let r2 = ((x as f64 + 0.5 - cx).powi(2) + (y as f64 + 0.5 - cy).powi(2)) / (m * m);
+                        let g = 1.0 + r2 * (k[0] + r2 * (k[1] + r2 * (k[2] + r2 * (k[3] + r2 * k[4]))));
+                        for p in 0..cpp {
+                            buf[(y * w + x) * cpp + p] *= g as f32;
+                        }
+                    }
+                }
+            }
+            Opcode::FixBadPixelsConstant { constant, .. } => {
+                let c = (*constant as f64 / 65535.0 * scale as f64) as f32;
+                let bad: Vec<(usize, usize)> =
+                    (0..h).flat_map(|y| (0..w).map(move |x| (x, y))).filter(|&(x, y)| cpp == 1 && buf[y * w + x] == c).collect();
+                fix_pixels(buf, w, h, cpp, cfa, &bad);
+            }
+            Opcode::FixBadPixelsList { points, rects, .. } => {
+                let mut bad: Vec<(usize, usize)> = points.iter().map(|&(r, c)| (c as usize, r as usize)).filter(|&(x, y)| x < w && y < h).collect();
+                for r in rects {
+                    let [t, l, b, rr] = r.map(|v| v as usize);
+                    for y in t..b.min(h) {
+                        for x in l..rr.min(w) {
+                            bad.push((x, y));
+                        }
+                    }
+                    if bad.len() > 1 << 22 {
+                        break;
+                    }
+                }
+                fix_pixels(buf, w, h, cpp, cfa, &bad);
+            }
+            _ => {}
+        }
+    }
+}
+
+/// Replace listed pixels with the mean of the nearest same-colour neighbours that are not themselves bad.
+fn fix_pixels(buf: &mut [f32], w: usize, h: usize, cpp: usize, cfa: Option<&Cfa>, bad: &[(usize, usize)]) {
+    if bad.is_empty() {
+        return;
+    }
+    let mut is_bad = std::collections::HashSet::with_capacity(bad.len());
+    for &b in bad {
+        is_bad.insert(b);
+    }
+    let step: isize = match cfa {
+        Some(c) if c.width == 2 && c.height == 2 => 2,
+        Some(c) => c.width.max(c.height) as isize,
+        None => 1,
+    };
+    for &(x, y) in bad {
+        for p in 0..cpp {
+            let mut sum = 0.0;
+            let mut n = 0;
+            for (dx, dy) in [(-step, 0), (step, 0), (0, -step), (0, step), (-step, -step), (step, step), (-step, step), (step, -step)] {
+                let (nx, ny) = (x as isize + dx, y as isize + dy);
+                if nx < 0 || ny < 0 || nx >= w as isize || ny >= h as isize || is_bad.contains(&(nx as usize, ny as usize)) {
+                    continue;
+                }
+                sum += buf[(ny as usize * w + nx as usize) * cpp + p];
+                n += 1;
+            }
+            if n > 0 {
+                buf[(y * w + x) * cpp + p] = sum / n as f32;
+            }
+        }
+    }
+}
+
+/// Apply `OpcodeList3` (on the demosaiced RGB image, values normalised to [0, 1]).
+pub fn apply_list3(list: &[Opcode], img: &mut Rgb32f) {
+    if list.is_empty() {
+        return;
+    }
+    let (w, h) = (img.width, img.height);
+    for op in list {
+        match op {
+            Opcode::WarpRectilinear { planes, center } => *img = warp_rectilinear(img, planes, *center),
+            Opcode::WarpFisheye { .. } | Opcode::TrimBounds { .. } | Opcode::Unknown { .. } => {}
+            other => {
+                let mut flat: Vec<f32> = img.data.iter().flat_map(|p| *p).collect();
+                apply_list(std::slice::from_ref(other), &mut flat, w, h, 3, None, 1.0);
+                for (d, c) in img.data.iter_mut().zip(flat.chunks_exact(3)) {
+                    *d = [c[0], c[1], c[2]];
+                }
+            }
+        }
+    }
+}
+
+/// DNG `WarpRectilinear`: for each output pixel, the source position is the radial + tangential model around the
+/// optical centre, with coordinates normalised by the distance to the farthest corner.
+pub fn warp_rectilinear(img: &Rgb32f, planes: &[[f64; 6]], center: [f64; 2]) -> Rgb32f {
+    let (w, h) = (img.width, img.height);
+    let (cx, cy) = (center[0] * (w as f64 - 1.0), center[1] * (h as f64 - 1.0));
+    let m = max_radius(cx, cy, w as f64 - 1.0, h as f64 - 1.0);
+    let mut out = Rgb32f::new(w, h);
+    let coef = |p: usize| planes[p.min(planes.len() - 1)];
+    lightcraft_raster::par_rows(&mut out.data, w, |y, row| {
+        for (x, px) in row.iter_mut().enumerate() {
+            for (c, v) in px.iter_mut().enumerate() {
+                let [kr0, kr1, kr2, kr3, kt0, kt1] = coef(c);
+                let dx = (x as f64 - cx) / m;
+                let dy = (y as f64 - cy) / m;
+                let r2 = dx * dx + dy * dy;
+                let f = kr0 + r2 * (kr1 + r2 * (kr2 + r2 * kr3));
+                let sx = dx * f + kt0 * 2.0 * dx * dy + kt1 * (r2 + 2.0 * dx * dx);
+                let sy = dy * f + kt1 * 2.0 * dx * dy + kt0 * (r2 + 2.0 * dy * dy);
+                let (fx, fy) = (cx + m * sx, cy + m * sy);
+                *v = img.sample_bilinear(fx as f32 + 0.5, fy as f32 + 0.5)[c];
+            }
+        }
+    });
+    out
+}
+
+#[cfg(test)]
+pub(crate) mod tests {
+    use super::*;
+
+    pub(crate) struct ListWriter(pub Vec<u8>, pub u32);
+
+    impl ListWriter {
+        pub fn new() -> Self {
+            ListWriter(Vec::new(), 0)
+        }
+        pub fn op(&mut self, id: u32, params: &[u8]) -> &mut Self {
+            for v in [id, 0x01030000, 0, params.len() as u32] {
+                self.0.extend_from_slice(&v.to_be_bytes());
+            }
+            self.0.extend_from_slice(params);
+            self.1 += 1;
+            self
+        }
+        pub fn finish(&self) -> Vec<u8> {
+            let mut v = self.1.to_be_bytes().to_vec();
+            v.extend_from_slice(&self.0);
+            v
+        }
+    }
+
+    pub(crate) fn be(vals: &[f64]) -> Vec<u8> {
+        vals.iter().flat_map(|v| v.to_be_bytes()).collect()
+    }
+    pub(crate) fn area(a: [u32; 8]) -> Vec<u8> {
+        a.iter().flat_map(|v| v.to_be_bytes()).collect()
+    }
+
+    #[test]
+    fn parse_and_apply_area_ops() {
+        let (w, h) = (4usize, 3usize);
+        let mut gm = area([0, 0, 3, 4, 0, 1, 1, 1]);
+        for v in [2u32, 2] {
+            gm.extend_from_slice(&v.to_be_bytes());
+        }
+        gm.extend(be(&[1.0, 1.0, 0.0, 0.0]));
+        gm.extend_from_slice(&1u32.to_be_bytes());
+        for g in [1.0f32, 2.0, 1.0, 2.0] {
+            gm.extend_from_slice(&g.to_be_bytes());
+        }
+        let mut dpr = area([0, 0, 3, 4, 0, 1, 1, 1]);
+        dpr.extend_from_slice(&3u32.to_be_bytes());
+        for d in [0.0f32, 0.1, 0.2] {
+            dpr.extend_from_slice(&d.to_be_bytes());
+        }
+        let mut poly = area([0, 0, 3, 4, 0, 1, 1, 2]);
+        poly.extend_from_slice(&1u32.to_be_bytes());
+        poly.extend(be(&[0.0, 2.0]));
+        let blob = ListWriter::new().op(9, &gm).op(10, &dpr).op(8, &poly).op(99, &[1, 2, 3]).finish();
+        let list = parse_list(&blob);
+        assert_eq!(list.len(), 4);
+        assert!(matches!(list[3], Opcode::Unknown { id: 99, .. }));
+        assert!(list[..3].iter().all(|o| o.is_applied()));
+        let mut buf = vec![0.25f32; w * h];
+        apply_list(&list, &mut buf, w, h, 1, None, 1.0);
+        // gain map: horizontal ramp 1 → 2 across relative x 0..1 (x / w)
+        let expect = |x: usize, y: usize| {
+            let g = 1.0 + x as f32 / w as f32;
+            let v = 0.25 * g + [0.0, 0.1, 0.2][y];
+            if x.is_multiple_of(2) { v * 2.0 } else { v }
+        };
+        for y in 0..h {
+            for x in 0..w {
+                assert!((buf[y * w + x] - expect(x, y)).abs() < 1e-6, "{x},{y}: {} vs {}", buf[y * w + x], expect(x, y));
+            }
+        }
+    }
+
+    #[test]
+    fn map_table_raw_scale_and_bad_pixels() {
+        let mut t = area([0, 0, 2, 2, 0, 1, 1, 1]);
+        t.extend_from_slice(&3u32.to_be_bytes());
+        for v in [0u16, 30000, 65535] {
+            t.extend_from_slice(&v.to_be_bytes());
+        }
+        let list = parse_list(&ListWriter::new().op(7, &t).finish());
+        let mut buf = vec![0.0, 1.0 / 65535.0, 2.0 / 65535.0, 0.5];
+        apply_list(&list, &mut buf, 2, 2, 1, None, 1.0);
+        assert_eq!(buf[0], 0.0);
+        assert!((buf[1] - 30000.0 / 65535.0).abs() < 1e-7);
+        assert_eq!(buf[3], 1.0);
+
+        let mut p = vec![];
+        for v in [7u32, 0] {
+            p.extend_from_slice(&v.to_be_bytes());
+        }
+        let list = parse_list(&ListWriter::new().op(4, &p).finish());
+        let cfa = Cfa::bayer("RGGB").unwrap();
+        let mut img: Vec<f32> = (0..36).map(|i| 100.0 + (i % 6) as f32).collect();
+        img[2 * 6 + 2] = 7.0;
+        apply_list(&list, &mut img, 6, 6, 1, Some(&cfa), 65535.0);
+        assert_eq!(img[2 * 6 + 2], 102.0);
+    }
+
+    #[test]
+    fn vignette_and_warp_identity() {
+        let mut img = Rgb32f::from_fn(9, 7, |x, y| [x as f32 / 9.0, y as f32 / 7.0, 0.5]);
+        let orig = img.clone();
+        let mut p = be(&[0.0; 5]);
+        p.extend(be(&[0.5, 0.5]));
+        let mut warp = 1u32.to_be_bytes().to_vec();
+        warp.extend(be(&[1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5]));
+        let list = parse_list(&ListWriter::new().op(3, &p).op(1, &warp).finish());
+        assert_eq!(list.len(), 2);
+        apply_list3(&list, &mut img);
+        for (a, b) in img.data.iter().zip(&orig.data) {
+            for c in 0..3 {
+                assert!((a[c] - b[c]).abs() < 1e-5);
+            }
+        }
+        // positive k0 brightens the corners more than the centre
+        let mut p = be(&[1.0, 0.0, 0.0, 0.0, 0.0]);
+        p.extend(be(&[0.5, 0.5]));
+        let list = parse_list(&ListWriter::new().op(3, &p).finish());
+        let mut flat = Rgb32f::filled(9, 9, [0.5; 3]);
+        apply_list3(&list, &mut flat);
+        assert!(flat.get(0, 0)[0] > 0.85 && (flat.get(4, 4)[0] - 0.5).abs() < 0.01, "{:?}", flat.get(0, 0));
+    }
+
+    #[test]
+    fn warp_barrel_moves_pixels_outward() {
+        let img = Rgb32f::from_fn(31, 31, |x, y| if x == 24 && y == 15 { [1.0; 3] } else { [0.0; 3] });
+        // source radius = r (1 + k r²) with k < 0 → content moves outward
+        let out = warp_rectilinear(&img, &[[1.0, -0.6, 0.0, 0.0, 0.0, 0.0]], [0.5, 0.5]);
+        let (bx, _) = (0..31).map(|x| (x, out.get(x, 15)[0])).fold((0, 0.0), |a, b| if b.1 > a.1 { b } else { a });
+        assert!(bx >= 25, "{bx}");
+    }
+
+    #[test]
+    fn write_parse_roundtrip() {
+        let a = Area { top: 1, left: 2, bottom: 30, right: 40, plane: 0, planes: 3, row_pitch: 1, col_pitch: 2 };
+        let list = vec![
+            Opcode::WarpRectilinear { planes: vec![[1.0, 0.01, -0.002, 0.0, 1e-4, -2e-4]; 3], center: [0.5, 0.49] },
+            Opcode::WarpFisheye { planes: vec![[1.0, 0.1, 0.0, 0.0]], center: [0.5, 0.5] },
+            Opcode::FixVignetteRadial { k: [0.1, 0.2, 0.0, 0.0, 0.0], center: [0.5, 0.5] },
+            Opcode::FixBadPixelsConstant { constant: 0, bayer_phase: 1 },
+            Opcode::FixBadPixelsList { bayer_phase: 0, points: vec![(3, 4)], rects: vec![[1, 1, 2, 2]] },
+            Opcode::TrimBounds { top: 1, left: 1, bottom: 9, right: 9 },
+            Opcode::MapTable { area: a, table: vec![0, 100, 65535] },
+            Opcode::MapPolynomial { area: a, coefficients: vec![0.0, 1.0, -0.1] },
+            Opcode::GainMap { area: a, points_v: 2, points_h: 1, spacing: [1.0, 1.0], origin: [0.0, 0.0], map_planes: 1, gains: vec![1.0, 1.5] },
+            Opcode::DeltaPerRow { area: a, deltas: vec![0.1, 0.2] },
+            Opcode::DeltaPerColumn { area: a, deltas: vec![0.3] },
+            Opcode::ScalePerRow { area: a, scales: vec![1.1] },
+            Opcode::ScalePerColumn { area: a, scales: vec![0.9, 1.0] },
+            Opcode::Unknown { id: 77, flags: 1, params: vec![1, 2, 3] },
+        ];
+        assert_eq!(parse_list(&write_list(&list)), list);
+    }
+
+    #[test]
+    fn truncated_lists_do_not_panic() {
+        let mut gm = area([0, 0, 3, 4, 0, 1, 1, 1]);
+        gm.extend_from_slice(&[0xff; 64]);
+        let blob = ListWriter::new().op(9, &gm).op(5, &[0xff; 12]).op(1, &[0, 0, 0, 9]).finish();
+        for n in 0..blob.len() {
+            let l = parse_list(&blob[..n]);
+            let mut buf = vec![0.5f32; 12];
+            apply_list(&l, &mut buf, 4, 3, 1, None, 1.0);
+        }
+    }
+}
diff --git a/crates/raw/src/preview.rs b/crates/raw/src/preview.rs
new file mode 100644
index 00000000..8b7aa21e
--- /dev/null
+++ b/crates/raw/src/preview.rs
@@ -0,0 +1,121 @@
+//! Embedded preview extraction: the largest baseline/progressive JPEG stored in a TIFF-based raw (IFD strips
+//! with JPEG compression, `JPEGInterchangeFormat` pointers in any IFD, Nikon/others' maker-note preview IFDs).
+
+use lightcraft_tiff::image::chunk_bytes;
+use lightcraft_tiff::tags as t;
+use lightcraft_tiff::{Ifd, Tiff, makernote};
+
+/// Whether `b` looks like a displayable (DCT) JPEG: SOI, and the first SOF marker is not lossless.
+fn is_dct_jpeg(b: &[u8]) -> bool {
+    if b.len() < 4 || b[0] != 0xff || b[1] != 0xd8 {
+        return false;
+    }
+    let mut i = 2;
+    while i + 4 <= b.len() {
+        if b[i] != 0xff {
+            return false;
+        }
+        let m = b[i + 1];
+        if m == 0xff {
+            i += 1;
+            continue;
+        }
+        match m {
+            0xc0..=0xc2 => return true,
+            0xc3 | 0xc5..=0xc7 | 0xcb | 0xcd..=0xcf => return false,
+            0xda | 0xd9 => return false,
+            _ => {}
+        }
+        let len = u16::from_be_bytes([b[i + 2], b[i + 3]]) as usize;
+        i += 2 + len;
+    }
+    false
+}
+
+fn candidates<'a>(data: &'a [u8], ifd: &Ifd, base: u64, out: &mut Vec<&'a [u8]>) {
+    if let (Some(off), Some(len)) = (ifd.u64(t::JPEG_INTERCHANGE_FORMAT), ifd.u64(t::JPEG_INTERCHANGE_FORMAT_LENGTH)) {
+        let off = off.saturating_add(base);
+        if let Some(s) = data.get(off as usize..(off.saturating_add(len) as usize).min(data.len())) {
+            out.push(s);
+        }
+    }
+    if matches!(ifd.u16(t::COMPRESSION), Some(6) | Some(7) | Some(34892))
+        && let Ok(info) = ifd.image()
+    {
+        let chunks = info.chunks(data.len() as u64);
+        if chunks.len() == 1
+            && let Some(s) = chunk_bytes(data, &chunks[0])
+        {
+            out.push(s);
+        }
+    }
+}
+
+/// The largest embedded JPEG preview, if any.
+pub fn embedded_preview(bytes: &[u8]) -> Option<Vec<u8>> {
+    let tiff = Tiff::parse(bytes).ok()?;
+    let mut found: Vec<&[u8]> = Vec::new();
+    for ifd in tiff.all_ifds() {
+        candidates(bytes, ifd, 0, &mut found);
+    }
+    // maker-note preview IFDs (e.g. Nikon PreviewIFD 0x0011 holds JPEGInterchangeFormat relative to the note base)
+    if let Some(exif) = tiff.exif()
+        && let Some(e) = exif.get(t::MAKER_NOTE)
+    {
+        let make = tiff.find(t::MAKE).and_then(|e| e.value.as_str()).unwrap_or_default();
+        if let Some(mn) = makernote::parse_makernote(bytes, e.offset, e.count() as u64, tiff.order, &make) {
+            candidates(bytes, &mn.ifd, mn.base, &mut found);
+            if let Some(off) = mn.ifd.u64(0x0011)
+                && let Ok((pifd, _)) = lightcraft_tiff::parse_ifd_at(bytes, mn.base + off, mn.order, mn.base, false, &Default::default())
+            {
+                candidates(bytes, &pifd, mn.base, &mut found);
+            }
+        }
+    }
+    found.into_iter().filter(|s| is_dct_jpeg(s)).max_by_key(|s| s.len()).map(|s| {
+        // trim trailing garbage after EOI when the length over-reports
+        let end = s.windows(2).rposition(|w| w == [0xff, 0xd9]).map(|p| p + 2).unwrap_or(s.len());
+        s[..end].to_vec()
+    })
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use lightcraft_tiff::{IfdBuilder, ImageData, TiffWriter, Value};
+
+    fn fake_jpeg(n: usize) -> Vec<u8> {
+        let mut j = vec![0xff, 0xd8, 0xff, 0xc0, 0x00, 0x0b, 8, 0, 1, 0, 1, 1, 1, 0x11, 0];
+        j.extend(std::iter::repeat_n(0x55u8, n));
+        j.extend_from_slice(&[0xff, 0xd9]);
+        j
+    }
+
+    #[test]
+    fn picks_largest_dct_jpeg() {
+        let small = fake_jpeg(10);
+        let big = fake_jpeg(500);
+        let lossless = crate::ljpeg::encode(&[1u16; 64], 8, 8, 1, 12, 1, 0);
+        let mut ifd0 = IfdBuilder::new();
+        ifd0.set(t::COMPRESSION, Value::Short(vec![6]));
+        ifd0.set(t::IMAGE_WIDTH, Value::Long(vec![1]));
+        ifd0.set(t::IMAGE_LENGTH, Value::Long(vec![1]));
+        ifd0.set_image(ImageData::Strips { rows_per_strip: 1, strips: vec![small.clone()] });
+        let mut raw = IfdBuilder::new();
+        raw.set(t::COMPRESSION, Value::Short(vec![7]));
+        raw.set(t::IMAGE_WIDTH, Value::Long(vec![8]));
+        raw.set(t::IMAGE_LENGTH, Value::Long(vec![8]));
+        raw.set_image(ImageData::Strips { rows_per_strip: 8, strips: vec![lossless] });
+        ifd0.add_sub_ifd(raw);
+        let mut sub = IfdBuilder::new();
+        sub.set(t::COMPRESSION, Value::Short(vec![7]));
+        sub.set(t::IMAGE_WIDTH, Value::Long(vec![2]));
+        sub.set(t::IMAGE_LENGTH, Value::Long(vec![2]));
+        sub.set(t::NEW_SUBFILE_TYPE, Value::Long(vec![1]));
+        sub.set_image(ImageData::Strips { rows_per_strip: 2, strips: vec![big.clone()] });
+        ifd0.add_sub_ifd(sub);
+        let bytes = TiffWriter::default().write(&[ifd0]).unwrap();
+        assert_eq!(embedded_preview(&bytes).unwrap(), big);
+        assert!(embedded_preview(b"nope").is_none());
+    }
+}
diff --git a/crates/raw/src/tiffraw.rs b/crates/raw/src/tiffraw.rs
new file mode 100644
index 00000000..fbf7e922
--- /dev/null
+++ b/crates/raw/src/tiffraw.rs
@@ -0,0 +1,226 @@
+//! Reading the pixel data of a TIFF IFD (strips or tiles; uncompressed, lossless JPEG, Deflate), in parallel.
+
+use crate::unpack::*;
+use crate::{MAX_SAMPLES, RawData, RawError, Result, ljpeg};
+use lightcraft_tiff::image::{Chunk, ImageInfo, chunk_bytes};
+use lightcraft_tiff::{ByteOrder, tags::compression as comp};
+use rayon::prelude::*;
+
+/// Bit packing for uncompressed integer data with bit depths other than 8/16.
+#[derive(Clone, Copy, Debug, PartialEq, Eq)]
+#[allow(dead_code)] // vendor decoders use the other variants
+pub enum Packing {
+    /// TIFF standard: MSB-first, rows start on byte boundaries.
+    Msb,
+    /// LSB-first little-endian bit stream (some vendor raws), rows byte aligned.
+    Lsb,
+    /// Each sample stored in 16 bits in file byte order, regardless of `bits`.
+    Word16,
+}
+
+enum ChunkPx {
+    U16(Vec<u16>),
+    F32(Vec<f32>),
+}
+
+/// Decode all chunks of `info` into one buffer of `width × height × cpp` samples.
+pub fn read_image(data: &[u8], info: &ImageInfo, order: ByteOrder, packing: Packing) -> Result<RawData> {
+    let (w, h) = (info.width as usize, info.height as usize);
+    let cpp = info.samples_per_pixel as usize;
+    let total = w.checked_mul(h).and_then(|v| v.checked_mul(cpp)).ok_or(RawError::Limit("image too large"))?;
+    if total > MAX_SAMPLES || total == 0 {
+        return Err(RawError::Limit("image too large"));
+    }
+    let bits = info.bits() as u32;
+    let float = info.sample_format == 3;
+    if float && !matches!(bits, 16 | 24 | 32) {
+        return Err(RawError::Unsupported(format!("{bits}-bit floating point samples")));
+    }
+    if !float && !(1..=16).contains(&bits) {
+        return Err(RawError::Unsupported(format!("{bits}-bit integer samples")));
+    }
+    let chunks = info.chunks(data.len() as u64);
+    if chunks.is_empty() {
+        return Err(RawError::Corrupt("no image data chunks".into()));
+    }
+    // Plausibility bound before allocating: no supported coding stores more than ~2000 samples per byte
+    // (Deflate of constant data is the extreme), so tiny files cannot trigger huge allocations.
+    let available: u64 = chunks.iter().map(|c| chunk_bytes(data, c).map_or(0, |s| s.len() as u64)).sum();
+    if available.saturating_mul(2048) < total as u64 {
+        return Err(RawError::Corrupt(format!("{available} bytes of image data cannot hold {total} samples")));
+    }
+    let planar = info.planar == 2 && cpp > 1;
+    let ccpp = if planar { 1 } else { cpp };
+    let decoded: Vec<Result<(Chunk, ChunkPx)>> =
+        chunks.par_iter().map(|c| decode_chunk(data, info, c, order, packing, ccpp, bits, float).map(|p| (*c, p))).collect();
+    let mut out_u16 = if float { Vec::new() } else { vec![0u16; total] };
+    let mut out_f32 = if float { vec![0f32; total] } else { Vec::new() };
+    let mut ok = 0usize;
+    let mut first_err = None;
+    for r in decoded {
+        let (c, px) = match r {
+            Ok(v) => v,
+            Err(e) => {
+                first_err.get_or_insert(e);
+                continue;
+            }
+        };
+        ok += 1;
+        let (cw, ch) = (c.width as usize, c.height as usize);
+        let (x0, y0) = (c.x as usize, c.y as usize);
+        if x0 >= w || y0 >= h {
+            continue;
+        }
+        let copy_w = cw.min(w - x0);
+        let copy_h = ch.min(h - y0);
+        for y in 0..copy_h {
+            for x in 0..copy_w {
+                for s in 0..ccpp {
+                    let src = (y * cw + x) * ccpp + s;
+                    let dst_s = if planar { c.plane as usize } else { s };
+                    if dst_s >= cpp {
+                        continue;
+                    }
+                    let dst = ((y0 + y) * w + x0 + x) * cpp + dst_s;
+                    match &px {
+                        ChunkPx::U16(v) => {
+                            if let (Some(o), Some(&s)) = (out_u16.get_mut(dst), v.get(src)) {
+                                *o = s;
+                            }
+                        }
+                        ChunkPx::F32(v) => {
+                            if let (Some(o), Some(&s)) = (out_f32.get_mut(dst), v.get(src)) {
+                                *o = s;
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+    if ok == 0 {
+        return Err(first_err.unwrap_or_else(|| RawError::Corrupt("no decodable chunks".into())));
+    }
+    Ok(if float { RawData::F32(out_f32) } else { RawData::U16(out_u16) })
+}
+
+#[allow(clippy::too_many_arguments)]
+fn decode_chunk(data: &[u8], info: &ImageInfo, c: &Chunk, order: ByteOrder, packing: Packing, cpp: usize, bits: u32, float: bool) -> Result<ChunkPx> {
+    let src = chunk_bytes(data, c).ok_or_else(|| RawError::Corrupt("chunk offset past end of file".into()))?;
+    let (cw, ch) = (c.width as usize, c.height as usize);
+    let n = cw.checked_mul(ch).and_then(|v| v.checked_mul(cpp)).ok_or(RawError::Limit("chunk too large"))?;
+    if n > MAX_SAMPLES {
+        return Err(RawError::Limit("chunk too large"));
+    }
+    match info.compression {
+        comp::NONE => unpack_chunk(src, order, packing, cw, ch, cpp, bits, float, 1),
+        comp::JPEG => {
+            let f = ljpeg::decode(src, n.saturating_mul(2).max(1 << 16))?;
+            if f.data.len() < n {
+                // Some writers encode edge tiles smaller than nominal: accept a frame that covers the rows it has.
+                if f.data.is_empty() || f.data.len() % (cw * cpp) != 0 {
+                    return Err(RawError::Corrupt(format!("lossless JPEG tile has {} samples, expected {n}", f.data.len())));
+                }
+            }
+            let mut v = f.data;
+            v.resize(n, 0);
+            Ok(ChunkPx::U16(v))
+        }
+        comp::ADOBE_DEFLATE | comp::DEFLATE => {
+            let bytes_per = if float { bits.div_ceil(8) as usize } else { 0 };
+            let row_bytes = if float { cw * cpp * bytes_per } else { (cw * cpp * bits as usize).div_ceil(8) };
+            let raw = inflate(src, row_bytes * ch)?;
+            unpack_chunk(&raw, order, packing, cw, ch, cpp, bits, float, info.predictor)
+        }
+        other => Err(RawError::Unsupported(format!("TIFF compression {other}"))),
+    }
+}
+
+#[allow(clippy::too_many_arguments)]
+fn unpack_chunk(
+    src: &[u8],
+    order: ByteOrder,
+    packing: Packing,
+    cw: usize,
+    ch: usize,
+    cpp: usize,
+    bits: u32,
+    float: bool,
+    predictor: u16,
+) -> Result<ChunkPx> {
+    let row_n = cw * cpp;
+    let (factor, fp) = match predictor {
+        1 => (0, false),
+        2 => (1, false),
+        3 => (1, true),
+        34892 => (2, false),
+        34893 => (4, false),
+        34894 => (2, true),
+        34895 => (4, true),
+        p => return Err(RawError::Unsupported(format!("predictor {p}"))),
+    };
+    if float {
+        let bp = bits.div_ceil(8) as usize;
+        let row_bytes = row_n * bp;
+        let mut out = vec![0f32; row_n * ch];
+        let mut rowbuf = vec![0u8; row_bytes];
+        for y in 0..ch {
+            let s = src.get(y * row_bytes..).unwrap_or(&[]);
+            let len = s.len().min(row_bytes);
+            rowbuf[..len].copy_from_slice(&s[..len]);
+            rowbuf[len..].fill(0);
+            let big = if fp {
+                undo_float_predictor(&mut rowbuf, row_n, bp, cpp * factor);
+                true
+            } else {
+                if factor > 0 {
+                    return Err(RawError::Unsupported("integer predictor on float data".into()));
+                }
+                order == ByteOrder::Big
+            };
+            for (i, o) in out[y * row_n..(y + 1) * row_n].iter_mut().enumerate() {
+                let b = &rowbuf[i * bp..(i + 1) * bp];
+                *o = match (bp, big) {
+                    (2, true) => f16_to_f32(u16::from_be_bytes([b[0], b[1]])),
+                    (2, false) => f16_to_f32(u16::from_le_bytes([b[0], b[1]])),
+                    (3, true) => f24_to_f32(u32::from_be_bytes([0, b[0], b[1], b[2]])),
+                    (3, false) => f24_to_f32(u32::from_le_bytes([b[0], b[1], b[2], 0])),
+                    (4, true) => f32::from_be_bytes([b[0], b[1], b[2], b[3]]),
+                    _ => f32::from_le_bytes([b[0], b[1], b[2], b[3]]),
+                };
+            }
+        }
+        return Ok(ChunkPx::F32(out));
+    }
+    if fp {
+        return Err(RawError::Unsupported("floating-point predictor on integer data".into()));
+    }
+    let mut out = vec![0u16; row_n * ch];
+    let (row_bytes, kind) = match (bits, packing) {
+        (8, _) => (row_n, 0),
+        (16, _) | (_, Packing::Word16) => (row_n * 2, 1),
+        (_, Packing::Msb) => ((row_n * bits as usize).div_ceil(8), 2),
+        (_, Packing::Lsb) => ((row_n * bits as usize).div_ceil(8), 3),
+    };
+    for y in 0..ch {
+        let s = src.get(y * row_bytes..).unwrap_or(&[]);
+        let s = &s[..s.len().min(row_bytes)];
+        let row = &mut out[y * row_n..(y + 1) * row_n];
+        match kind {
+            0 => row.iter_mut().zip(s).for_each(|(o, &b)| *o = b as u16),
+            1 => read_u16s(s, order, row),
+            2 => unpack_msb(s, bits, row),
+            _ => unpack_lsb(s, bits, row),
+        }
+        if factor > 0 {
+            if bits == 8 {
+                for i in cpp * factor..row.len() {
+                    row[i] = (row[i] as u8).wrapping_add(row[i - cpp * factor] as u8) as u16;
+                }
+            } else {
+                undo_diff_u16(row, cpp * factor);
+            }
+        }
+    }
+    Ok(ChunkPx::U16(out))
+}
diff --git a/crates/raw/src/unpack.rs b/crates/raw/src/unpack.rs
new file mode 100644
index 00000000..6d86e02f
--- /dev/null
+++ b/crates/raw/src/unpack.rs
@@ -0,0 +1,186 @@
+//! Sample unpacking: bit packing, byte orders, 16/24/32-bit floats, Deflate and TIFF / DNG predictors.
+//!
+//! Sources: TIFF 6.0 §14 (Differencing Predictor), Adobe Photoshop TIFF Technical Note 3 (floating-point
+//! predictor), DNG 1.7 (predictors 34892–34895, bit packing, 16/24-bit floats), RFC 1950/1951 (zlib/deflate).
+
+use crate::{RawError, Result};
+use lightcraft_tiff::ByteOrder;
+
+/// Unpack `out.len()` MSB-first packed samples of `bits` (1..=16) from `src` (TIFF bit order).
+pub fn unpack_msb(src: &[u8], bits: u32, out: &mut [u16]) {
+    debug_assert!((1..=16).contains(&bits));
+    let mut acc: u64 = 0;
+    let mut nacc: u32 = 0;
+    let mut it = src.iter();
+    let mask = (1u64 << bits) - 1;
+    for o in out.iter_mut() {
+        while nacc < bits {
+            acc = (acc << 8) | *it.next().unwrap_or(&0) as u64;
+            nacc += 8;
+        }
+        *o = ((acc >> (nacc - bits)) & mask) as u16;
+        nacc -= bits;
+    }
+}
+
+/// Unpack LSB-first packed samples (bits taken from the low end of little-endian words).
+pub fn unpack_lsb(src: &[u8], bits: u32, out: &mut [u16]) {
+    let mut acc: u64 = 0;
+    let mut nacc: u32 = 0;
+    let mut it = src.iter();
+    let mask = (1u64 << bits) - 1;
+    for o in out.iter_mut() {
+        while nacc < bits {
+            acc |= (*it.next().unwrap_or(&0) as u64) << nacc;
+            nacc += 8;
+        }
+        *o = (acc & mask) as u16;
+        acc >>= bits;
+        nacc -= bits;
+    }
+}
+
+pub fn read_u16s(src: &[u8], order: ByteOrder, out: &mut [u16]) {
+    for (o, c) in out.iter_mut().zip(src.chunks_exact(2)) {
+        *o = order.u16([c[0], c[1]]);
+    }
+}
+
+/// IEEE 754 half → f32.
+pub fn f16_to_f32(h: u16) -> f32 {
+    let sign = ((h >> 15) as u32) << 31;
+    let exp = ((h >> 10) & 0x1f) as u32;
+    let man = (h & 0x3ff) as u32;
+    let bits = if exp == 0 {
+        if man == 0 {
+            sign
+        } else {
+            // subnormal: normalise
+            let mut e = 127 - 15 + 1;
+            let mut m = man;
+            while m & 0x400 == 0 {
+                m <<= 1;
+                e -= 1;
+            }
+            sign | (e << 23) | ((m & 0x3ff) << 13)
+        }
+    } else if exp == 31 {
+        sign | 0x7f80_0000 | (man << 13)
+    } else {
+        sign | ((exp + 127 - 15) << 23) | (man << 13)
+    };
+    f32::from_bits(bits)
+}
+
+/// DNG 24-bit float (1 sign, 7 exponent bias 63, 16 mantissa) → f32.
+pub fn f24_to_f32(v: u32) -> f32 {
+    let sign = (v >> 23) & 1;
+    let exp = (v >> 16) & 0x7f;
+    let man = v & 0xffff;
+    let bits = if exp == 0 && man == 0 {
+        sign << 31
+    } else if exp == 0x7f {
+        (sign << 31) | 0x7f80_0000 | (man << 7)
+    } else if exp == 0 {
+        // denormal: value = man * 2^(1-63-16)
+        let v = man as f32 * 2f32.powi(1 - 63 - 16);
+        return if sign == 1 { -v } else { v };
+    } else {
+        (sign << 31) | ((exp + 127 - 63) << 23) | (man << 7)
+    };
+    f32::from_bits(bits)
+}
+
+/// Inflate a zlib stream, refusing to produce more than `limit` bytes.
+pub fn inflate(src: &[u8], limit: usize) -> Result<Vec<u8>> {
+    miniz_oxide::inflate::decompress_to_vec_zlib_with_limit(src, limit).map_err(|e| RawError::Corrupt(format!("deflate: {e:?}")))
+}
+
+/// Undo horizontal differencing on one row of integer samples with a sample `stride` (cpp × factor).
+pub fn undo_diff_u16(row: &mut [u16], stride: usize) {
+    for i in stride..row.len() {
+        row[i] = row[i].wrapping_add(row[i - stride]);
+    }
+}
+
+pub fn undo_diff_u8(row: &mut [u8], stride: usize) {
+    for i in stride..row.len() {
+        row[i] = row[i].wrapping_add(row[i - stride]);
+    }
+}
+
+/// Undo the floating-point predictor on one row: `row` holds `n` samples of `bytes_per` bytes, byte-differenced
+/// with stride `stride` samples, then stored as byte planes (most significant first). Returns the samples as
+/// big-endian byte groups in place.
+pub fn undo_float_predictor(row: &mut [u8], n: usize, bytes_per: usize, stride: usize) {
+    undo_diff_u8(row, stride);
+    let tmp = row.to_vec();
+    for i in 0..n {
+        for b in 0..bytes_per {
+            if let (Some(dst), Some(&s)) = (row.get_mut(i * bytes_per + b), tmp.get(b * n + i)) {
+                *dst = s;
+            }
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn packing() {
+        // 12-bit MSB: 0xABC, 0xDEF -> AB CD EF
+        let mut o = [0u16; 2];
+        unpack_msb(&[0xab, 0xcd, 0xef], 12, &mut o);
+        assert_eq!(o, [0xabc, 0xdef]);
+        unpack_lsb(&[0xbc, 0xfa, 0xde], 12, &mut o);
+        assert_eq!(o, [0xabc, 0xdef]);
+        let mut o = [0u16; 3];
+        unpack_msb(&[0b1010_1100, 0b0100_0000], 3, &mut o);
+        assert_eq!(o, [0b101, 0b011, 0b000]);
+        let mut o = [0u16; 2];
+        read_u16s(&[1, 2, 3, 4], ByteOrder::Big, &mut o);
+        assert_eq!(o, [0x0102, 0x0304]);
+        // short source: pads with zeros
+        let mut o = [9u16; 4];
+        unpack_msb(&[0xff], 14, &mut o);
+        assert_eq!(o, [0x3fc0, 0, 0, 0]);
+    }
+
+    #[test]
+    fn floats() {
+        assert_eq!(f16_to_f32(0x3c00), 1.0);
+        assert_eq!(f16_to_f32(0xc000), -2.0);
+        assert_eq!(f16_to_f32(0x0001), 2f32.powi(-24));
+        assert!(f16_to_f32(0x7c00).is_infinite());
+        // 24-bit: 1.0 = exp 63, mantissa 0
+        assert_eq!(f24_to_f32(63 << 16), 1.0);
+        assert_eq!(f24_to_f32((1 << 23) | (64 << 16) | 0x8000), -3.0);
+        assert_eq!(f24_to_f32(0), 0.0);
+    }
+
+    #[test]
+    fn predictors() {
+        let mut r = [1u16, 1, 1, 1, 65535, 2];
+        undo_diff_u16(&mut r, 2);
+        assert_eq!(r, [1, 1, 2, 2, 1, 4]);
+        // float predictor: two 32-bit samples 1.0 (3f800000), 2.0 (40000000), stride 1
+        let planes: [u8; 8] = [0x3f, 0x40, 0x80, 0x00, 0x00, 0x00, 0x00, 0x00];
+        let mut row = planes;
+        for i in (1..row.len()).rev() {
+            row[i] = row[i].wrapping_sub(row[i - 1]);
+        }
+        undo_float_predictor(&mut row, 2, 4, 1);
+        assert_eq!(f32::from_be_bytes([row[0], row[1], row[2], row[3]]), 1.0);
+        assert_eq!(f32::from_be_bytes([row[4], row[5], row[6], row[7]]), 2.0);
+    }
+
+    #[test]
+    fn inflate_limit() {
+        let z = miniz_oxide::deflate::compress_to_vec_zlib(&[7u8; 1000], 6);
+        assert_eq!(inflate(&z, 1000).unwrap(), vec![7u8; 1000]);
+        assert!(inflate(&z, 10).is_err());
+        assert!(inflate(&[1, 2, 3], 10).is_err());
+    }
+}
diff --git a/crates/raw/src/vendor/arw.rs b/crates/raw/src/vendor/arw.rs
new file mode 100644
index 00000000..b6197a79
--- /dev/null
+++ b/crates/raw/src/vendor/arw.rs
@@ -0,0 +1,5 @@
+use crate::{RawError, RawImage, Result};
+
+pub(crate) fn decode(_bytes: &[u8]) -> Result<RawImage> {
+    Err(RawError::Unsupported("arw: not yet".into()))
+}
diff --git a/crates/raw/src/vendor/cr2.rs b/crates/raw/src/vendor/cr2.rs
new file mode 100644
index 00000000..7452793d
--- /dev/null
+++ b/crates/raw/src/vendor/cr2.rs
@@ -0,0 +1,5 @@
+use crate::{RawError, RawImage, Result};
+
+pub(crate) fn decode(_bytes: &[u8]) -> Result<RawImage> {
+    Err(RawError::Unsupported("cr2: not yet".into()))
+}
diff --git a/crates/raw/src/vendor/mod.rs b/crates/raw/src/vendor/mod.rs
new file mode 100644
index 00000000..54cea8b6
--- /dev/null
+++ b/crates/raw/src/vendor/mod.rs
@@ -0,0 +1,5 @@
+//! Vendor raw formats (TIFF-based).
+
+pub mod arw;
+pub mod cr2;
+pub mod nef;
diff --git a/crates/raw/src/vendor/nef.rs b/crates/raw/src/vendor/nef.rs
new file mode 100644
index 00000000..b82f3453
--- /dev/null
+++ b/crates/raw/src/vendor/nef.rs
@@ -0,0 +1,5 @@
+use crate::{RawError, RawImage, Result};
+
+pub(crate) fn decode(_bytes: &[u8]) -> Result<RawImage> {
+    Err(RawError::Unsupported("nef: not yet".into()))
+}
diff --git a/crates/raw/tests/dng.rs b/crates/raw/tests/dng.rs
new file mode 100644
index 00000000..6fa7b670
--- /dev/null
+++ b/crates/raw/tests/dng.rs
@@ -0,0 +1,426 @@
+//! Synthetic DNG round trips: every storage variant we decode is written with our own TIFF writer (and LJ92
+//! encoder) and must decode bit-exactly.
+
+use lightcraft_color::{D65, Mat3, Xy};
+use lightcraft_raw::opcodes::{Area, Opcode};
+use lightcraft_raw::{
+    BlackLevel, Cfa, ColorData, DngCompression, DngWriteOptions, Method, OpcodeLists, Orientation, RawData, RawFormat, RawImage, Rect, decode,
+    embedded_preview, probe, write_dng,
+};
+use lightcraft_tiff::tags::{self as t, compression, photometric};
+use lightcraft_tiff::{ByteOrder, IfdBuilder, ImageData, TiffWriter, Value};
+
+fn scene_value(x: usize, y: usize, c: u8, max: f32) -> f32 {
+    let base = [0.55, 0.8, 0.4][c as usize];
+    let v = base * (0.2 + 0.7 * ((x as f32 * 0.13).sin() * 0.5 + 0.5) * ((y as f32 * 0.07).cos() * 0.3 + 0.7));
+    v * max
+}
+
+fn synthetic(w: usize, h: usize, cfa: Option<Cfa>, cpp: usize) -> RawImage {
+    let black = 256.0;
+    let white = 16000.0;
+    let mut data = Vec::with_capacity(w * h * cpp);
+    for y in 0..h {
+        for x in 0..w {
+            for s in 0..cpp {
+                let c = match &cfa {
+                    Some(p) => p.color_at(x, y),
+                    None => s as u8,
+                };
+                data.push((black + scene_value(x, y, c, white - black)).round() as u16);
+            }
+        }
+    }
+    let color = ColorData {
+        illuminant: [17, 21],
+        color_matrix: [
+            Some(Mat3([[0.9, 0.2, -0.15], [-0.3, 1.25, 0.08], [0.02, -0.12, 0.85]])),
+            Some(Mat3([[0.7, 0.3, -0.1], [-0.35, 1.3, 0.1], [0.05, -0.2, 1.0]])),
+        ],
+        as_shot_neutral: Some([0.5, 1.0, 0.7]),
+        baseline_exposure: 0.35,
+        ..Default::default()
+    };
+    let metadata = lightcraft_meta::Metadata { make: Some("Synth".into()), model: Some("Cam 1".into()), rating: Some(3), ..Default::default() };
+    RawImage {
+        format: RawFormat::Dng,
+        width: w,
+        height: h,
+        cpp,
+        data: RawData::U16(data),
+        cfa,
+        bits: 14,
+        black: BlackLevel::uniform(black),
+        white: vec![white],
+        active_area: Rect::new(2, 1, w - 4, h - 2),
+        crop: Rect::new(3, 2, w - 12, h - 8),
+        orientation: Orientation::Rotate90,
+        color,
+        wb_multipliers: None,
+        linearized: false,
+        opcodes: OpcodeLists::default(),
+        metadata,
+    }
+}
+
+fn assert_same(a: &RawImage, b: &RawImage) {
+    assert_eq!((a.width, a.height, a.cpp), (b.width, b.height, b.cpp));
+    assert_eq!(a.data, b.data);
+    assert_eq!(a.cfa, b.cfa);
+    assert_eq!(a.active_area, b.active_area);
+    assert_eq!(a.crop, b.crop);
+    assert_eq!(a.orientation, b.orientation);
+    assert!((0..a.cpp).all(|s| a.white_at(s) == b.white_at(s)));
+    assert_eq!(a.black.values, b.black.values.iter().take(a.black.values.len()).copied().collect::<Vec<_>>());
+    for i in 0..2 {
+        let (x, y) = (a.color.color_matrix[i].unwrap(), b.color.color_matrix[i].unwrap());
+        for r in 0..3 {
+            for c in 0..3 {
+                assert!((x.0[r][c] - y.0[r][c]).abs() < 1e-6);
+            }
+        }
+    }
+    assert_eq!(a.color.illuminant, b.color.illuminant);
+    let (n1, n2) = (a.color.as_shot_neutral.unwrap(), b.color.as_shot_neutral.unwrap());
+    assert!((0..3).all(|i| (n1[i] - n2[i]).abs() < 1e-6));
+    assert!((a.color.baseline_exposure - b.color.baseline_exposure).abs() < 1e-6);
+    assert_eq!(b.metadata.make.as_deref(), Some("Synth"));
+    assert_eq!(b.metadata.rating, Some(3));
+}
+
+#[test]
+fn writer_roundtrip_all_layouts() {
+    for (w, h) in [(64, 48), (70, 45), (33, 17)] {
+        for cfa in [Some(Cfa::bayer("RGGB").unwrap()), Some(Cfa::bayer("GBRG").unwrap()), Some(Cfa::xtrans()), None] {
+            let cpp = if cfa.is_some() { 1 } else { 3 };
+            let raw = synthetic(w, h, cfa, cpp);
+            for comp in [DngCompression::Uncompressed, DngCompression::Lj92 { tile: 16 }, DngCompression::Lj92 { tile: 256 }] {
+                for order in [ByteOrder::Little, ByteOrder::Big] {
+                    let bytes = write_dng(&raw, &DngWriteOptions { compression: comp, order, xmp: None }).unwrap();
+                    assert_eq!(probe(&bytes), Some(RawFormat::Dng));
+                    let back = decode(&bytes).unwrap();
+                    assert_same(&raw, &back);
+                }
+            }
+        }
+    }
+}
+
+#[test]
+fn normalized_and_develop() {
+    let mut raw = synthetic(40, 30, Some(Cfa::bayer("RGGB").unwrap()), 1);
+    let n = raw.normalized().unwrap();
+    assert_eq!((n.width, n.height), (36, 28));
+    // the CFA is re-anchored at the active area origin (2, 1): RGGB shifted by (2,1) → GBRG
+    assert_eq!(n.cfa.as_ref().unwrap().name(), "GBRG");
+    let RawData::U16(d) = &raw.data else { unreachable!() };
+    let expect = (d[40 + 2] as f32 - 256.0) / (16000.0 - 256.0);
+    assert!((n.data[0] - expect).abs() < 1e-6);
+    let rgb = raw.develop(Method::Ahd).unwrap();
+    assert_eq!((rgb.width, rgb.height), (28, 22));
+    // constant scene → constant output, exactly 0.5 after normalisation
+    if let RawData::U16(d) = &mut raw.data {
+        d.iter_mut().for_each(|v| *v = 256 + (16000 - 256) / 2);
+    }
+    for m in [Method::Bilinear, Method::Ppg, Method::Ahd] {
+        let rgb = raw.develop(m).unwrap();
+        assert!(rgb.data.iter().all(|p| p.iter().all(|v| (v - 0.5).abs() < 1e-4)), "{m:?}");
+    }
+}
+
+// ---------------------------------------------------------------- hand-built variants
+
+fn base_ifd(w: usize, h: usize, bits: u16, cpp: usize, cfa: bool) -> IfdBuilder {
+    let mut ifd = IfdBuilder::new();
+    ifd.set(t::NEW_SUBFILE_TYPE, Value::Long(vec![0]));
+    ifd.set(t::IMAGE_WIDTH, Value::Long(vec![w as u32]));
+    ifd.set(t::IMAGE_LENGTH, Value::Long(vec![h as u32]));
+    ifd.set(t::BITS_PER_SAMPLE, Value::Short(vec![bits; cpp]));
+    ifd.set(t::SAMPLES_PER_PIXEL, Value::Short(vec![cpp as u16]));
+    ifd.set(t::PHOTOMETRIC, Value::Short(vec![if cfa { photometric::CFA } else { photometric::LINEAR_RAW }]));
+    if cfa {
+        ifd.set(t::CFA_REPEAT_PATTERN_DIM, Value::Short(vec![2, 2]));
+        ifd.set(t::CFA_PATTERN_EP, Value::Byte(vec![1, 0, 2, 1]));
+    }
+    ifd
+}
+
+fn dng(raw: IfdBuilder, order: ByteOrder) -> Vec<u8> {
+    let mut ifd0 = IfdBuilder::new();
+    ifd0.set(t::NEW_SUBFILE_TYPE, Value::Long(vec![1]));
+    ifd0.set(t::DNG_VERSION, Value::Byte(vec![1, 6, 0, 0]));
+    ifd0.set(t::MAKE, Value::Ascii("Hand".into()));
+    ifd0.set(t::COLOR_MATRIX_1, Value::SRational(vec![(1, 1), (0, 1), (0, 1), (0, 1), (1, 1), (0, 1), (0, 1), (0, 1), (1, 1)]));
+    ifd0.set(t::CALIBRATION_ILLUMINANT_1, Value::Short(vec![21]));
+    ifd0.set(t::AS_SHOT_WHITE_XY, Value::Rational(vec![(3127, 10000), (3290, 10000)]));
+    // an 8-bit RGB thumbnail in IFD0 (must not be mistaken for the raw)
+    ifd0.set(t::IMAGE_WIDTH, Value::Long(vec![4]));
+    ifd0.set(t::IMAGE_LENGTH, Value::Long(vec![2]));
+    ifd0.set(t::BITS_PER_SAMPLE, Value::Short(vec![8, 8, 8]));
+    ifd0.set(t::SAMPLES_PER_PIXEL, Value::Short(vec![3]));
+    ifd0.set(t::PHOTOMETRIC, Value::Short(vec![photometric::RGB]));
+    ifd0.set_image(ImageData::Strips { rows_per_strip: 2, strips: vec![vec![128; 24]] });
+    let mut exif = IfdBuilder::new();
+    exif.set(t::ISO_SPEED, Value::Short(vec![200]));
+    ifd0.set_child(t::EXIF_IFD, exif);
+    ifd0.add_sub_ifd(raw);
+    TiffWriter::new(order, false).write(&[ifd0]).unwrap()
+}
+
+fn pattern(w: usize, h: usize, cpp: usize, bits: u32) -> Vec<u16> {
+    let max = (1u32 << bits) - 1;
+    (0..w * h * cpp).map(|i| ((i as u32).wrapping_mul(2654435761u32) >> 7) % (max + 1)).map(|v| v as u16).collect()
+}
+
+fn pack_msb(v: &[u16], bits: u32, w: usize) -> Vec<u8> {
+    let mut out = Vec::new();
+    for row in v.chunks(w) {
+        let mut acc = 0u64;
+        let mut n = 0;
+        for &s in row {
+            acc = (acc << bits) | s as u64;
+            n += bits;
+            while n >= 8 {
+                out.push((acc >> (n - 8)) as u8);
+                n -= 8;
+            }
+        }
+        if n > 0 {
+            out.push((acc << (8 - n)) as u8);
+        }
+    }
+    out
+}
+
+fn u16_bytes(v: &[u16], order: ByteOrder) -> Vec<u8> {
+    let mut b = Vec::new();
+    v.iter().for_each(|&x| order.put_u16(&mut b, x));
+    b
+}
+
+#[test]
+fn packed_bit_depths() {
+    let (w, h) = (13, 9);
+    for bits in [10u32, 12, 14] {
+        let px = pattern(w, h, 1, bits);
+        let mut raw = base_ifd(w, h, bits as u16, 1, true);
+        raw.set(t::COMPRESSION, Value::Short(vec![1]));
+        let strips: Vec<Vec<u8>> = px.chunks(w * 4).map(|c| pack_msb(c, bits, w)).collect();
+        raw.set_image(ImageData::Strips { rows_per_strip: 4, strips });
+        let img = decode(&dng(raw, ByteOrder::Big)).unwrap();
+        assert_eq!(img.data, RawData::U16(px), "{bits} bits");
+        assert_eq!(img.white, vec![((1u32 << bits) - 1) as f32]);
+        assert_eq!(img.cfa.unwrap().name(), "GRBG");
+        assert_eq!(img.metadata.iso, Some(200));
+        assert_eq!(img.color.as_shot_white_xy, Some(Xy::new(0.3127, 0.329)));
+    }
+}
+
+#[test]
+fn deflate_with_integer_predictors() {
+    let (w, h, cpp) = (10, 6, 3);
+    let px = pattern(w, h, cpp, 16);
+    for (pred, factor) in [(1u16, 0usize), (2, 1), (34892, 2), (34893, 4)] {
+        for order in [ByteOrder::Little, ByteOrder::Big] {
+            let mut raw = base_ifd(w, h, 16, cpp, false);
+            raw.set(t::COMPRESSION, Value::Short(vec![compression::ADOBE_DEFLATE]));
+            raw.set(t::PREDICTOR, Value::Short(vec![pred]));
+            raw.set_image(ImageData::Tiles {
+                tile_width: 16,
+                tile_height: 16,
+                tiles: vec![{
+                    // one 16×16 tile: re-lay the 10×6 image into the tile grid
+                    let mut tile = vec![0u16; 16 * 16 * cpp];
+                    for y in 0..h {
+                        for x in 0..w * cpp {
+                            tile[y * 16 * cpp + x] = px[y * w * cpp + x];
+                        }
+                    }
+                    let mut te = tile.clone();
+                    if factor > 0 {
+                        for row in te.chunks_mut(16 * cpp) {
+                            for i in (cpp * factor..row.len()).rev() {
+                                row[i] = row[i].wrapping_sub(row[i - cpp * factor]);
+                            }
+                        }
+                    }
+                    miniz_oxide::deflate::compress_to_vec_zlib(&u16_bytes(&te, order), 6)
+                }],
+            });
+            let img = decode(&dng(raw, order)).unwrap();
+            assert_eq!(img.data, RawData::U16(px.clone()), "predictor {pred} {order:?}");
+            assert!(img.cfa.is_none());
+            assert_eq!(img.cpp, 3);
+        }
+    }
+}
+
+fn fp_encode(vals: &[f32], bytes_per: usize, stride: usize, w_samples: usize) -> Vec<u8> {
+    // big-endian byte planes per row, then byte differencing
+    let mut out = Vec::new();
+    for row in vals.chunks(w_samples) {
+        let be: Vec<[u8; 4]> = row
+            .iter()
+            .map(|&v| match bytes_per {
+                4 => v.to_be_bytes(),
+                2 => {
+                    let h = f32_to_f16(v).to_be_bytes();
+                    [h[0], h[1], 0, 0]
+                }
+                _ => {
+                    let b = f32_to_f24(v).to_be_bytes();
+                    [b[1], b[2], b[3], 0]
+                }
+            })
+            .collect();
+        let mut planes = vec![0u8; row.len() * bytes_per];
+        for (i, b) in be.iter().enumerate() {
+            for k in 0..bytes_per {
+                planes[k * row.len() + i] = b[k];
+            }
+        }
+        for i in (stride..planes.len()).rev() {
+            planes[i] = planes[i].wrapping_sub(planes[i - stride]);
+        }
+        out.extend(planes);
+    }
+    out
+}
+
+fn f32_to_f16(v: f32) -> u16 {
+    // exact for the test values (multiples of 1/64 in [0, 4))
+    let b = v.to_bits();
+    let sign = ((b >> 16) & 0x8000) as u16;
+    if v == 0.0 {
+        return sign;
+    }
+    let exp = ((b >> 23) & 0xff) as i32 - 127 + 15;
+    let man = ((b >> 13) & 0x3ff) as u16;
+    sign | ((exp as u16) << 10) | man
+}
+
+fn f32_to_f24(v: f32) -> u32 {
+    let b = v.to_bits();
+    if v == 0.0 {
+        return 0;
+    }
+    let sign = (b >> 31) << 23;
+    let exp = (((b >> 23) & 0xff) as i32 - 127 + 63) as u32;
+    sign | (exp << 16) | ((b >> 7) & 0xffff)
+}
+
+#[test]
+fn floating_point_dngs() {
+    let (w, h) = (12, 5);
+    let vals: Vec<f32> = (0..w * h).map(|i| ((i * 37) % 200) as f32 / 64.0).collect();
+    for (bits, pred, stride_factor) in [(32u16, 3u16, 1usize), (32, 34894, 2), (32, 34895, 4), (16, 3, 1), (24, 3, 1), (24, 34894, 2)] {
+        let bp = bits as usize / 8;
+        let mut raw = base_ifd(w, h, bits, 1, true);
+        raw.set(t::SAMPLE_FORMAT, Value::Short(vec![3]));
+        raw.set(t::COMPRESSION, Value::Short(vec![compression::DEFLATE]));
+        raw.set(t::PREDICTOR, Value::Short(vec![pred]));
+        let z = miniz_oxide::deflate::compress_to_vec_zlib(&fp_encode(&vals, bp, stride_factor, w), 6);
+        raw.set_image(ImageData::Strips { rows_per_strip: h as u32, strips: vec![z] });
+        let img = decode(&dng(raw, ByteOrder::Little)).unwrap();
+        assert_eq!(img.data, RawData::F32(vals.clone()), "{bits}-bit predictor {pred}");
+        assert_eq!(img.white, vec![1.0]);
+    }
+    // uncompressed half floats, big-endian file
+    let mut raw = base_ifd(w, h, 16, 1, true);
+    raw.set(t::SAMPLE_FORMAT, Value::Short(vec![3]));
+    raw.set(t::COMPRESSION, Value::Short(vec![1]));
+    let bytes: Vec<u8> = vals.iter().flat_map(|&v| f32_to_f16(v).to_be_bytes()).collect();
+    raw.set_image(ImageData::Strips { rows_per_strip: h as u32, strips: vec![bytes] });
+    assert_eq!(decode(&dng(raw, ByteOrder::Big)).unwrap().data, RawData::F32(vals));
+}
+
+#[test]
+fn lj92_strips_and_levels_and_linearization() {
+    let (w, h) = (20, 12);
+    let px = pattern(w, h, 1, 12);
+    let mut raw = base_ifd(w, h, 12, 1, true);
+    raw.set(t::COMPRESSION, Value::Short(vec![compression::JPEG]));
+    let strips: Vec<Vec<u8>> = px.chunks(w * 5).map(|c| lightcraft_raw::ljpeg::encode(c, w, c.len() / w, 1, 12, 6, 0)).collect();
+    raw.set_image(ImageData::Strips { rows_per_strip: 5, strips });
+    let table: Vec<u16> = (0..4096u32).map(|i| (i * i / 4096) as u16).collect();
+    raw.set(t::LINEARIZATION_TABLE, Value::Short(table.clone()));
+    raw.set(t::BLACK_LEVEL_REPEAT_DIM, Value::Short(vec![2, 2]));
+    raw.set(t::BLACK_LEVEL, Value::Rational(vec![(10, 1), (11, 1), (12, 1), (13, 1)]));
+    raw.set(t::BLACK_LEVEL_DELTA_H, Value::SRational(vec![(1, 2); 16]));
+    raw.set(t::BLACK_LEVEL_DELTA_V, Value::SRational(vec![(-1, 1); 10]));
+    raw.set(t::WHITE_LEVEL, Value::Short(vec![4000]));
+    raw.set(t::ACTIVE_AREA, Value::Long(vec![1, 2, 11, 18]));
+    raw.set(t::DEFAULT_CROP_ORIGIN, Value::Rational(vec![(1, 1), (2, 1)]));
+    raw.set(t::DEFAULT_CROP_SIZE, Value::Rational(vec![(12, 1), (6, 1)]));
+    let img = decode(&dng(raw, ByteOrder::Little)).unwrap();
+    let expect: Vec<u16> = px.iter().map(|&v| table[v as usize]).collect();
+    assert_eq!(img.data, RawData::U16(expect.clone()));
+    assert!(img.linearized);
+    assert_eq!(img.active_area, Rect::new(2, 1, 16, 10));
+    assert_eq!(img.crop, Rect::new(1, 2, 12, 6));
+    assert_eq!(img.black.at(1, 1, 0, 1), 13.0 + 0.5 - 1.0);
+    let n = img.normalized().unwrap();
+    let v = expect[(1 + 1) * w + 2 + 1] as f32;
+    let b = 13.0 + 0.5 - 1.0;
+    assert!((n.data[16 + 1] - (v - b) / (4000.0 - 11.5)).abs() < 1e-5);
+    let rgb = img.develop(Method::Ppg).unwrap();
+    assert_eq!((rgb.width, rgb.height), (12, 6));
+}
+
+#[test]
+fn opcodes_are_parsed_and_applied() {
+    let mut raw = synthetic(32, 24, Some(Cfa::bayer("RGGB").unwrap()), 1);
+    let area = Area { top: 0, left: 0, bottom: 24, right: 32, plane: 0, planes: 1, row_pitch: 1, col_pitch: 1 };
+    raw.opcodes.list2 = vec![Opcode::ScalePerRow { area, scales: vec![2.0; 24] }];
+    raw.opcodes.list3 = vec![Opcode::FixVignetteRadial { k: [0.0; 5], center: [0.5, 0.5] }, Opcode::Unknown { id: 99, flags: 1, params: vec![] }];
+    let bytes = write_dng(&raw, &DngWriteOptions::default()).unwrap();
+    let back = decode(&bytes).unwrap();
+    assert_eq!(back.opcodes, raw.opcodes);
+    let plain = {
+        let mut r = raw.clone();
+        r.opcodes = Default::default();
+        r.normalized().unwrap()
+    };
+    let scaled = back.normalized().unwrap();
+    for (a, b) in plain.data.iter().zip(&scaled.data) {
+        assert!((b - 2.0 * a).abs() < 1e-5);
+    }
+}
+
+#[test]
+fn preview_and_probe() {
+    let raw = synthetic(16, 16, Some(Cfa::bayer("RGGB").unwrap()), 1);
+    let bytes = write_dng(&raw, &DngWriteOptions::default()).unwrap();
+    // the raw tiles are lossless JPEG: never offered as a preview
+    assert!(embedded_preview(&bytes).is_none());
+    assert_eq!(probe(&bytes), Some(RawFormat::Dng));
+}
+
+#[test]
+fn colour_transform_from_decoded_dng() {
+    let raw = decode(&write_dng(&synthetic(16, 16, Some(Cfa::bayer("RGGB").unwrap()), 1), &DngWriteOptions::default()).unwrap()).unwrap();
+    let xy = lightcraft_raw::color::as_shot_white_xy(&raw);
+    let t = lightcraft_raw::color::camera_transform(&raw, xy);
+    assert!(!t.matrix_is_fallback);
+    // as-shot neutral (0.5, 1, 0.7) → multipliers (2, 1, 1/0.7) up to normalisation
+    assert!((t.wb[0] / t.wb[1] - 2.0).abs() < 1e-3, "{:?}", t.wb);
+    assert!((t.wb[2] / t.wb[1] - 1.0 / 0.7).abs() < 1e-3);
+    assert!((t.baseline_exposure - 0.35).abs() < 1e-6);
+    let (m, _) = lightcraft_raw::color::camera_to_rec2020(&raw, D65);
+    let o = m.apply([1.0; 3]);
+    assert!(o.iter().all(|v| (v - 1.0).abs() < 1e-9));
+}
+
+#[test]
+fn rejects_unsupported_and_broken() {
+    let (w, h) = (8, 8);
+    let mut raw = base_ifd(w, h, 8, 1, true);
+    raw.set(t::COMPRESSION, Value::Short(vec![compression::LOSSY_JPEG]));
+    raw.set_image(ImageData::Strips { rows_per_strip: 8, strips: vec![vec![0; 10]] });
+    assert!(matches!(decode(&dng(raw, ByteOrder::Little)), Err(lightcraft_raw::RawError::Unsupported(_))));
+    let mut raw = base_ifd(w, h, 16, 1, true);
+    raw.set(t::CFA_PATTERN_EP, Value::Byte(vec![0, 1, 3, 1]));
+    raw.set(t::COMPRESSION, Value::Short(vec![1]));
+    raw.set_image(ImageData::Strips { rows_per_strip: 8, strips: vec![vec![0; 128]] });
+    assert!(decode(&dng(raw, ByteOrder::Little)).is_err());
+}
diff --git a/crates/raw/tests/robust.rs b/crates/raw/tests/robust.rs
new file mode 100644
index 00000000..7bd5a5e5
--- /dev/null
+++ b/crates/raw/tests/robust.rs
@@ -0,0 +1,104 @@
+//! Malformed-input robustness: decoding random or corrupted files must never panic.
+
+use lightcraft_raw::{
+    BlackLevel, Cfa, ColorData, DngCompression, DngWriteOptions, Method, OpcodeLists, Orientation, RawData, RawFormat, RawImage, Rect,
+};
+use proptest::prelude::*;
+
+fn sample(comp: DngCompression) -> Vec<u8> {
+    let (w, h) = (24usize, 16usize);
+    let data: Vec<u16> = (0..w * h).map(|i| (i * 97 % 4000) as u16).collect();
+    let raw = RawImage {
+        format: RawFormat::Dng,
+        width: w,
+        height: h,
+        cpp: 1,
+        data: RawData::U16(data),
+        cfa: Some(Cfa::bayer("RGGB").unwrap()),
+        bits: 12,
+        black: BlackLevel::uniform(64.0),
+        white: vec![4095.0],
+        active_area: Rect::new(0, 0, w, h),
+        crop: Rect::new(0, 0, w, h),
+        orientation: Orientation::Normal,
+        color: ColorData::default(),
+        wb_multipliers: None,
+        linearized: false,
+        opcodes: OpcodeLists::default(),
+        metadata: Default::default(),
+    };
+    lightcraft_raw::write_dng(&raw, &DngWriteOptions { compression: comp, ..Default::default() }).unwrap()
+}
+
+fn exercise(bytes: &[u8]) {
+    let _ = lightcraft_raw::probe(bytes);
+    let _ = lightcraft_raw::embedded_preview(bytes);
+    if let Ok(img) = lightcraft_raw::decode(bytes)
+        && img.width * img.height <= 1 << 16
+    {
+        let _ = img.develop(Method::Ahd);
+        let xy = lightcraft_raw::color::as_shot_white_xy(&img);
+        let _ = lightcraft_raw::color::camera_transform(&img, xy);
+    }
+}
+
+proptest! {
+    #![proptest_config(ProptestConfig { cases: 400, .. ProptestConfig::default() })]
+
+    #[test]
+    fn mutated_dngs_never_panic(flips in proptest::collection::vec((any::<usize>(), any::<u8>()), 1..24), cut in any::<usize>(), lj in any::<bool>()) {
+        let mut data = sample(if lj { DngCompression::Lj92 { tile: 16 } } else { DngCompression::Uncompressed });
+        for (i, v) in flips {
+            let n = data.len();
+            data[i % n] = v;
+        }
+        let keep = if cut % 3 == 0 { cut % data.len() } else { data.len() };
+        exercise(&data[..keep]);
+    }
+
+    #[test]
+    fn random_tiff_like_bytes_never_panic(mut data in proptest::collection::vec(any::<u8>(), 8..700), kind in 0usize..4) {
+        let heads: [&[u8]; 4] = [b"II*\0\x08\0\0\0", b"MM\0*\0\0\0\x08", b"II*\0\x10\0\0\0CR\x02\0", b"II+\0\x08\0\0\0"];
+        let h = heads[kind];
+        let n = h.len().min(data.len());
+        data[..n].copy_from_slice(&h[..n]);
+        exercise(&data);
+    }
+
+    #[test]
+    fn random_lj92_never_panics(data in proptest::collection::vec(any::<u8>(), 0..400)) {
+        let mut s = vec![0xff, 0xd8, 0xff, 0xc3, 0, 11, 12, 0, 8, 0, 8, 1, 1, 0x11, 0];
+        s.extend_from_slice(&data);
+        let _ = lightcraft_raw::ljpeg::decode(&s, 1 << 20);
+    }
+}
+
+#[test]
+fn truncation_at_every_length() {
+    for comp in [DngCompression::Uncompressed, DngCompression::Lj92 { tile: 16 }] {
+        let d = sample(comp);
+        for n in (0..d.len()).step_by(7) {
+            exercise(&d[..n]);
+        }
+    }
+}
+
+#[test]
+fn hostile_dimensions_are_rejected_quickly() {
+    // a DNG claiming a 60000 × 60000 16-bit image with 10 bytes of data
+    use lightcraft_tiff::tags as t;
+    use lightcraft_tiff::{IfdBuilder, ImageData, TiffWriter, Value};
+    let mut ifd = IfdBuilder::new();
+    ifd.set(t::DNG_VERSION, Value::Byte(vec![1, 4, 0, 0]));
+    ifd.set(t::IMAGE_WIDTH, Value::Long(vec![60000]));
+    ifd.set(t::IMAGE_LENGTH, Value::Long(vec![60000]));
+    ifd.set(t::BITS_PER_SAMPLE, Value::Short(vec![16]));
+    ifd.set(t::PHOTOMETRIC, Value::Short(vec![32803]));
+    ifd.set(t::CFA_PATTERN_EP, Value::Byte(vec![0, 1, 1, 2]));
+    ifd.set(t::COMPRESSION, Value::Short(vec![7]));
+    ifd.set_image(ImageData::Tiles { tile_width: 60000, tile_height: 60000, tiles: vec![vec![0xff, 0xd8, 0xff, 0xd9]] });
+    let bytes = TiffWriter::default().write(&[ifd]).unwrap();
+    let t0 = std::time::Instant::now();
+    assert!(lightcraft_raw::decode(&bytes).is_err());
+    assert!(t0.elapsed().as_secs_f64() < 2.0);
+}