shifathosn
/
lightcraft
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

M3.2-M3.6: lightcraft-raw DNG, LJ92, demosaic, highlights, DNG colour 

- LJ92 (T.81 process 14) decoder + Annex K encoder: predictors 1-7,
  1-4 components, 2-16 bit, point transform, restart intervals
- DNG: uncompressed (packed 1-16 bit, 16/24/32-bit float), LJ92 tiles/
  strips, Deflate with predictors 2/3/34892-34895; CFA + LinearRaw;
  linearization, BlackLevelRepeatDim/Delta, WhiteLevel, ActiveArea,
  DefaultCrop; colour tags; OpcodeList1-3 parse/write/apply
- DNG writer (uncompressed / LJ92 tiles) for export and round trips
- demosaic: bilinear, PPG, AHD (Hirakawa-Parks), X-Trans directional
- highlight clip/reconstruct; DNG colour model with dual illuminant,
  forward matrices, as-shot neutral -> xy, Rec.2020 transform
- embedded preview extraction; malformed-input proptests

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Brandon Thomas 2026-09-30 11:33:19 -04:00
parent cbe9d8989a
commit 05450c791c
23 changed files with 4715 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

27
Cargo.lock generated
View File

@ -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"

12
crates/raw/Cargo.toml
View File

@ -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 }

377
crates/raw/src/color.rs Normal file
View File

@ -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:?}");
}
}

165
crates/raw/src/demosaic/ahd.rs Normal file
View File

@ -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
}

78
crates/raw/src/demosaic/bilinear.rs Normal file
View File

@ -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
}

253
crates/raw/src/demosaic/mod.rs Normal file
View File

@ -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());
}
}

89
crates/raw/src/demosaic/ppg.rs Normal file
View File

@ -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
}

95
crates/raw/src/demosaic/xtrans.rs Normal file
View File

@ -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
}

167
crates/raw/src/dng.rs Normal file
View File

@ -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)
}

201
crates/raw/src/dngwrite.rs Normal file
View File

@ -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])?)
}

180
crates/raw/src/highlight.rs Normal file
View File

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

489
crates/raw/src/lib.rs
View File

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

786
crates/raw/src/ljpeg.rs Normal file
View File

@ -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());
}
}

713
crates/raw/src/opcodes.rs Normal file
View File

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

121
crates/raw/src/preview.rs Normal file
View File

@ -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());
}
}

226
crates/raw/src/tiffraw.rs Normal file
View File

@ -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))
}

186
crates/raw/src/unpack.rs Normal file
View File

@ -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());
}
}

5
crates/raw/src/vendor/arw.rs vendored Normal file
View File

@ -0,0 +1,5 @@
use crate::{RawError, RawImage, Result};
pub(crate) fn decode(_bytes: &[u8]) -> Result<RawImage> {
Err(RawError::Unsupported("arw: not yet".into()))
}

5
crates/raw/src/vendor/cr2.rs vendored Normal file
View File

@ -0,0 +1,5 @@
use crate::{RawError, RawImage, Result};
pub(crate) fn decode(_bytes: &[u8]) -> Result<RawImage> {
Err(RawError::Unsupported("cr2: not yet".into()))
}

5
crates/raw/src/vendor/mod.rs vendored Normal file
View File

@ -0,0 +1,5 @@
//! Vendor raw formats (TIFF-based).
pub mod arw;
pub mod cr2;
pub mod nef;

5
crates/raw/src/vendor/nef.rs vendored Normal file
View File

@ -0,0 +1,5 @@
use crate::{RawError, RawImage, Result};
pub(crate) fn decode(_bytes: &[u8]) -> Result<RawImage> {
Err(RawError::Unsupported("nef: not yet".into()))
}

426
crates/raw/tests/dng.rs Normal file
View File

@ -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());
}

104
crates/raw/tests/robust.rs Normal file
View File

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