fix(raw): decode Samsung linear DNGs correctly

rawler/src/decoders/dng.rs

@@ -18,6 +18,24 @@ pub struct DngDecoder<'a> {
   tiff: GenericTiffReader,
 }
 
+fn validate_blacklevel_count(level_count: usize, repeat: (usize, usize), cpp: usize) -> Result<()> {
+  let expected_levels = repeat
+    .0
+    .checked_mul(repeat.1)
+    .and_then(|count| count.checked_mul(cpp))
+    .ok_or_else(|| format!("BlackLevel repeat size overflow: {}x{} with {cpp} samples per pixel", repeat.0, repeat.1))?;
+  if level_count != expected_levels {
+    return Err(
+      format!(
+        "BlackLevel count mismatch: expected {expected_levels} values for repeat {}x{} with {cpp} samples per pixel, found {level_count}",
+        repeat.0, repeat.1
+      )
+      .into(),
+    );
+  }
+  Ok(())
+}
+
 impl<'a> DngDecoder<'a> {
   pub fn new(_file: &RawSource, tiff: GenericTiffReader, rawloader: &'a RawLoader) -> Result<DngDecoder<'a>> {
     Ok(DngDecoder { tiff, rawloader })
@@ -327,6 +345,7 @@ impl<'a> DngDecoder<'a> {
           log::warn!("File has BlackLevelRepeatDim tag but with invalid length: {}", value.len());
         }
       }
+      validate_blacklevel_count(levels.len(), repeat, cpp)?;
       Ok(Some(BlackLevel::new(&levels, repeat.1, repeat.0, cpp)))
     } else {
       Ok(None)
@@ -434,3 +453,19 @@ impl<'a> DngDecoder<'a> {
     Ok(result)
   }
 }
+
+#[cfg(test)]
+mod tests {
+  use super::validate_blacklevel_count;
+
+  #[test]
+  fn rejects_blacklevel_count_mismatch_before_construction() {
+    let error = validate_blacklevel_count(1, (2, 2), 3).unwrap_err();
+    assert!(error.to_string().contains("expected 12 values"));
+  }
+
+  #[test]
+  fn accepts_row_column_sample_blacklevel_count() {
+    validate_blacklevel_count(12, (2, 2), 3).unwrap();
+  }
+}

rawler/src/decompressors/ljpeg/decompressors.rs

@@ -21,19 +21,29 @@ pub fn decode_ljpeg(ljpeg: &LjpegDecompressor, out: &mut [u16], x: usize, stripw
   let htable = |index: usize| -> &HuffTable { &ljpeg.dhts[ljpeg.sof.components[index].dc_tbl_num] };
   let mut pump = BitPumpJPEG::new(ljpeg.buffer);
   let base_prediction = 1 << (ljpeg.sof.precision - ljpeg.point_transform - 1);
+  let total_mcus = ljpeg
+    .sof
+    .width
+    .checked_mul(height)
+    .ok_or_else(|| format!("ljpeg: MCU count overflow for {}x{height}", ljpeg.sof.width))?;
+  let mut next_restart_marker = 0;
 
   // initialize first pixel components
   for c in 0..ncomp {
     out[x + c] = (base_prediction + htable(c).huff_decode(&mut pump)?) as u16;
   }
+  let mut decoded_mcus = 1;
+  let mut reset_prediction = consume_restart_if_needed(ljpeg, &mut pump, decoded_mcus, total_mcus, &mut next_restart_marker)?;
 
   let skip_x = ljpeg.sof.width - width / ncomp;
 
   for row in 0..height {
     let startcol = if row == 0 { x + ncomp } else { x }; // skip first pixel in first row
     for col in (startcol..(width + x)).step_by(ncomp) {
       for c in 0..ncomp {
-        let p: i32 = if col == x {
+        let p: i32 = if reset_prediction {
+          base_prediction
+        } else if col == x {
           // At start of line predictor starts with start of previous line
           out[(row - 1) * stripwidth + x + c] as i32
         } else {
@@ -83,16 +93,38 @@ pub fn decode_ljpeg(ljpeg: &LjpegDecompressor, out: &mut [u16], x: usize, stripw
         let diff = htable(c).huff_decode(&mut pump)?;
         out[row * stripwidth + col + c] = (p + diff) as u16;
       }
+      decoded_mcus += 1;
+      reset_prediction = consume_restart_if_needed(ljpeg, &mut pump, decoded_mcus, total_mcus, &mut next_restart_marker)?;
     }
     for _ in 0..skip_x {
+      // This MCU is outside the requested output width, but it still counts
+      // towards the restart interval and its entropy data must be consumed.
       for c in 0..ncomp {
         // Skip extra encoded differences if the ljpeg frame is wider than the output
         htable(c).huff_decode(&mut pump)?;
       }
+      decoded_mcus += 1;
+      reset_prediction = consume_restart_if_needed(ljpeg, &mut pump, decoded_mcus, total_mcus, &mut next_restart_marker)?;
     }
   }
 
-  Ok(())
+  if height == ljpeg.sof.height { pump.validate_end_of_scan() } else { Ok(()) }
+}
+
+fn consume_restart_if_needed(
+  ljpeg: &LjpegDecompressor,
+  pump: &mut BitPumpJPEG<'_>,
+  decoded_mcus: usize,
+  total_mcus: usize,
+  next_restart_marker: &mut u8,
+) -> Result<bool, String> {
+  if ljpeg.restart_interval != 0 && decoded_mcus % ljpeg.restart_interval == 0 && decoded_mcus < total_mcus {
+    pump.consume_restart_marker(*next_restart_marker)?;
+    *next_restart_marker = (*next_restart_marker + 1) % 8;
+    Ok(true)
+  } else {
+    Ok(false)
+  }
 }
 
 fn set_yuv_420(out: &mut [u16], row: usize, col: usize, width: usize, y1: i32, y2: i32, y3: i32, y4: i32, cb: i32, cr: i32) {

rawler/src/decompressors/ljpeg/mod.rs

@@ -16,6 +16,7 @@ enum Marker {
   EOI = 0xd9,  // end of image
   SOS = 0xda,  // start of scan
   DQT = 0xdb,  // quantization tables
+  DRI = 0xdd,  // restart interval
   Fill = 0xff,
 }
 
@@ -128,6 +129,7 @@ pub struct LjpegDecompressor<'a> {
   sof: SOFInfo,
   predictor: usize,
   point_transform: usize,
+  restart_interval: usize,
   dhts: Vec<HuffTable>,
 }
 
@@ -148,6 +150,7 @@ impl<'a> LjpegDecompressor<'a> {
     let mut dht_huffval = [[0_u32; 256]; 4];
     let pred;
     let pt;
+    let mut restart_interval = 0;
     loop {
       let marker = LjpegDecompressor::get_next_marker(&mut input, true)?;
       if marker == m(Marker::SOF3) {
@@ -159,6 +162,8 @@ impl<'a> LjpegDecompressor<'a> {
       } else if marker == m(Marker::DHT) {
         // Huffman table settings
         LjpegDecompressor::parse_dht(&mut input, &mut dht_init, &mut dht_bits, &mut dht_huffval)?;
+      } else if marker == m(Marker::DRI) {
+        restart_interval = LjpegDecompressor::parse_dri(&mut input)?;
       } else if marker == m(Marker::SOS) {
         // Start of the actual stream, we can decode after this
         let (a, b) = sof.parse_sos(&mut input)?;
@@ -183,13 +188,14 @@ impl<'a> LjpegDecompressor<'a> {
     }
 
     log::debug!(
-      "LJPEGDecompressor: super_h: {}, super_v: {}, pred: {}, pt: {}, prec: {}, cps: {}",
+      "LJPEGDecompressor: super_h: {}, super_v: {}, pred: {}, pt: {}, prec: {}, cps: {}, restart_interval: {}",
       sof.components[0].super_h,
       sof.components[0].super_v,
       pred,
       pt,
       sof.precision,
       sof.cps,
+      restart_interval,
     );
 
     if sof.components[0].super_h == 2 && sof.components[0].super_v == 2 {
@@ -204,6 +210,7 @@ impl<'a> LjpegDecompressor<'a> {
       sof,
       predictor: pred,
       point_transform: pt,
+      restart_interval,
       dhts,
     })
   }
@@ -269,11 +276,23 @@ impl<'a> LjpegDecompressor<'a> {
     Ok(())
   }
 
+  fn parse_dri(input: &mut ByteStream) -> Result<usize, String> {
+    if input.remaining_bytes() < 4 {
+      return Err(format!("ljpeg: truncated DRI segment ({} bytes remain)", input.remaining_bytes()));
+    }
+    let length = input.get_u16();
+    if length != 4 {
+      return Err(format!("ljpeg: invalid DRI length {length}"));
+    }
+    Ok(input.get_u16() as usize)
+  }
+
   /// Handle special SONY YUV 4:2:0 encoding in ILCE-7RM5
   pub fn decode_sony(&self, out: &mut [u16], x: usize, stripwidth: usize, width: usize, height: usize, dummy: bool) -> Result<(), String> {
     if dummy {
       return Ok(());
     }
+    self.validate_restart_support()?;
     log::debug!("LJPEG decode with special Sony mode");
     if self.sof.components[0].super_h == 2 && self.sof.components[0].super_v == 2 {
       decode_sony_ljpeg_420(self, out, width, height)
@@ -297,6 +316,7 @@ impl<'a> LjpegDecompressor<'a> {
     if dummy {
       return Ok(());
     }
+    self.validate_restart_support()?;
 
     if self.sof.components[0].super_h == 2 && self.sof.components[0].super_v == 2 {
       decode_ljpeg_420(self, out, width, height)
@@ -349,6 +369,25 @@ impl<'a> LjpegDecompressor<'a> {
     )
   }
 
+  fn validate_restart_support(&self) -> Result<(), String> {
+    if self.restart_interval == 0 {
+      return Ok(());
+    }
+
+    if let Some(component) = self.sof.components.iter().find(|component| component.super_h != 1 || component.super_v != 1) {
+      return Err(format!(
+        "ljpeg: restart markers are not supported for component {} sampling {}x{}",
+        component.id, component.super_h, component.super_v
+      ));
+    }
+
+    if !(1..=7).contains(&self.predictor) {
+      return Err(format!("ljpeg: restart markers are not supported with predictor {}", self.predictor));
+    }
+
+    Ok(())
+  }
+
   pub fn width(&self) -> usize {
     self.sof.width * self.sof.cps
   }

rawler/src/imgop/raw.rs

@@ -2,9 +2,9 @@
 // Copyright 2021 Daniel Vogelbacher <daniel@chaospixel.com>
 
 use super::xyz::Illuminant;
-use crate::imgop::Rect;
 use crate::imgop::matrix::{multiply, normalize, pseudo_inverse};
 use crate::imgop::xyz::SRGB_TO_XYZ_D65;
+use crate::imgop::{Point, Rect};
 use crate::pixarray::{Color2D, RgbF32};
 use crate::rawimage::{BlackLevel, RawPhotometricInterpretation, WhiteLevel};
 
@@ -151,6 +151,140 @@ pub fn correct_blacklevel(raw: &mut [f32], blacklevel: &[f32], whitelevel: &[f32
   }
 }
 
+/// Correct black and white levels for an interleaved linear raw image.
+///
+/// DNG black levels are stored in row-column-sample order and may repeat over
+/// multiple pixels, while white levels are stored once per sample plane. The
+/// repeat pattern is relative to `blacklevel_origin` (the top-left corner of
+/// the DNG ActiveArea).
+pub fn correct_blacklevel_linear(
+  raw: &mut [f32],
+  width: usize,
+  height: usize,
+  cpp: usize,
+  blacklevel: &BlackLevel,
+  whitelevel: &WhiteLevel,
+  blacklevel_origin: Point,
+) -> crate::Result<()> {
+  if width == 0 || height == 0 || cpp == 0 {
+    return Err(format!("Invalid linear raw dimensions: {width}x{height} with {cpp} components per pixel").into());
+  }
+
+  let row_len = width
+    .checked_mul(cpp)
+    .ok_or_else(|| format!("Linear raw row size overflow: width {width}, cpp {cpp}"))?;
+  let expected_raw_len = row_len
+    .checked_mul(height)
+    .ok_or_else(|| format!("Linear raw image size overflow: {width}x{height} with {cpp} components per pixel"))?;
+  if raw.len() != expected_raw_len {
+    return Err(format!("Linear raw data length mismatch: expected {expected_raw_len}, found {}", raw.len()).into());
+  }
+
+  if blacklevel.width == 0 || blacklevel.height == 0 || blacklevel.cpp == 0 {
+    return Err(
+      format!(
+        "Invalid black level repeat dimensions: {}x{} with {} components per pixel",
+        blacklevel.width, blacklevel.height, blacklevel.cpp
+      )
+      .into(),
+    );
+  }
+  if blacklevel.cpp != 1 && blacklevel.cpp != cpp {
+    return Err(format!("Black level component count mismatch: expected 1 or {cpp}, found {}", blacklevel.cpp).into());
+  }
+
+  let expected_blacklevel_len = blacklevel
+    .width
+    .checked_mul(blacklevel.height)
+    .and_then(|count| count.checked_mul(blacklevel.cpp))
+    .ok_or_else(|| {
+      format!(
+        "Black level repeat size overflow: {}x{} with {} components per pixel",
+        blacklevel.width, blacklevel.height, blacklevel.cpp
+      )
+    })?;
+  if blacklevel.levels.len() != expected_blacklevel_len {
+    return Err(
+      format!(
+        "Black level data length mismatch: expected {expected_blacklevel_len}, found {}",
+        blacklevel.levels.len()
+      )
+      .into(),
+    );
+  }
+
+  let whitelevels = match whitelevel.0.as_slice() {
+    [level] => vec![*level as f32; cpp],
+    levels if levels.len() == cpp => levels.iter().map(|level| *level as f32).collect(),
+    levels => {
+      return Err(format!("White level component count mismatch: expected 1 or {cpp}, found {}", levels.len()).into());
+    }
+  };
+  let blacklevels = blacklevel.as_vec();
+  if blacklevels.iter().any(|level| !level.is_finite()) {
+    return Err("Black level contains a non-finite value".into());
+  }
+
+  // DNG normalization uses the maximum computed black level for each sample
+  // plane, even when the repeating pattern contains different local values.
+  let mut max_blacklevels = vec![f32::NEG_INFINITY; cpp];
+  for cell in blacklevels.chunks_exact(blacklevel.cpp) {
+    for channel in 0..cpp {
+      let black_channel = if blacklevel.cpp == 1 { 0 } else { channel };
+      max_blacklevels[channel] = max_blacklevels[channel].max(cell[black_channel]);
+    }
+  }
+  let scales: Vec<f32> = whitelevels
+    .iter()
+    .zip(&max_blacklevels)
+    .enumerate()
+    .map(|(channel, (white, black))| {
+      let scale = *white - *black;
+      if scale.is_finite() && scale > 0.0 {
+        Ok(scale)
+      } else {
+        Err(format!("Invalid black/white level range for channel {channel}: black {black}, white {white}"))
+      }
+    })
+    .collect::<std::result::Result<_, _>>()?;
+
+  // Most linear DNGs, including Samsung's 2x2 patterns, repeat identical
+  // values. Keep that common case on the vectorized per-channel path.
+  let first_cell = &blacklevels[..blacklevel.cpp];
+  if blacklevels.chunks_exact(blacklevel.cpp).all(|cell| cell == first_cell) {
+    let flat_blacklevels: Vec<f32> = (0..cpp).map(|channel| first_cell[if blacklevel.cpp == 1 { 0 } else { channel }]).collect();
+    correct_blacklevel(raw, &flat_blacklevels, &whitelevels);
+    return Ok(());
+  }
+
+  let origin_x = blacklevel_origin.x % blacklevel.width;
+  let origin_y = blacklevel_origin.y % blacklevel.height;
+  raw.par_chunks_exact_mut(row_len).enumerate().for_each(|(y, row)| {
+    let repeat_y = repeat_position(y, origin_y, blacklevel.height);
+    row.chunks_exact_mut(cpp).enumerate().for_each(|(x, pixel)| {
+      let repeat_x = repeat_position(x, origin_x, blacklevel.width);
+      let cell_offset = (repeat_y * blacklevel.width + repeat_x) * blacklevel.cpp;
+      for channel in 0..cpp {
+        let black_channel = if blacklevel.cpp == 1 { 0 } else { channel };
+        let corrected = pixel[channel] - blacklevels[cell_offset + black_channel];
+        pixel[channel] = if corrected.is_sign_negative() { 0.0 } else { corrected / scales[channel] };
+      }
+    });
+  });
+
+  Ok(())
+}
+
+#[inline]
+fn repeat_position(position: usize, origin: usize, repeat: usize) -> usize {
+  let position = position % repeat;
+  if position >= origin {
+    position - origin
+  } else {
+    repeat - (origin - position)
+  }
+}
+
 /// Correct data by blacklevel and whitelevel on CFA (bayer) data.
 ///
 /// The output is between 0.0 .. 1.0.