From 273907b3d1926f2d823c39de40ca21090cc1869c Mon Sep 17 00:00:00 2001 From: Marcin Zawalski Date: Sun, 25 Oct 2026 09:07:19 +0100 Subject: [PATCH 01/12] perf(core): EXIF orientation in Rust apply_exif_orientation turned full-resolution buffers with NumPy rot90/flip plus a contiguous copy on one thread. negpy_core.orient copies pixels of any byte size in parallel, in 64x64 tiles for the 90-degree turns. Output is the same bytes as the NumPy path for every orientation and dtype. 23.5 MP float32 RGB, 24 threads: orientation 6 328 -> 27 ms, orientation 3 171 -> 24 ms. Runs on frame open, export decode and twice in lens correction for rotated sources. --- .github/workflows/ci.yml | 2 +- negpy/kernel/image/logic.py | 21 ++----- negpy_core/src/lib.rs | 1 + negpy_core/src/orient.rs | 110 +++++++++++++++++++++++++++++++++++ negpy_core/src/py/convert.rs | 21 +++++++ tests/test_core_orient.py | 55 ++++++++++++++++++ 6 files changed, 193 insertions(+), 17 deletions(-) create mode 100644 negpy_core/src/orient.rs create mode 100644 tests/test_core_orient.py diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 3213b83d8..33f7cdc5a 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -129,4 +129,4 @@ jobs: run: uv sync --group dev --group plustek - name: negpy_core pins and oracles - run: uv run pytest tests/test_percentiles.py tests/test_core_same_bits.py tests/test_core_resize.py tests/test_core_filter.py tests/test_core_morph.py tests/test_core_gray.py tests/test_core_contours.py tests/test_core_shapes.py tests/test_core_components.py tests/test_core_color.py tests/test_core_warp.py tests/test_core_lab.py tests/test_core_altprocess.py + run: uv run pytest tests/test_percentiles.py tests/test_core_same_bits.py tests/test_core_orient.py tests/test_core_resize.py tests/test_core_filter.py tests/test_core_morph.py tests/test_core_gray.py tests/test_core_contours.py tests/test_core_shapes.py tests/test_core_components.py tests/test_core_color.py tests/test_core_warp.py tests/test_core_lab.py tests/test_core_altprocess.py diff --git a/negpy/kernel/image/logic.py b/negpy/kernel/image/logic.py index 9411d0970..109ef31d3 100644 --- a/negpy/kernel/image/logic.py +++ b/negpy/kernel/image/logic.py @@ -265,23 +265,12 @@ def apply_exif_orientation(arr: np.ndarray, orientation: Optional[int]) -> np.nd Bake an EXIF orientation value (1-8) into pixels so the array displays upright. Works on HxW (IR) and HxWxC (RGB) arrays. Returns the input unchanged for 1/None. """ - if not orientation or orientation == 1: + if orientation not in (2, 3, 4, 5, 6, 7, 8) or arr.size == 0: return arr - if orientation == 2: - return np.ascontiguousarray(np.fliplr(arr)) - if orientation == 3: - return np.ascontiguousarray(np.rot90(arr, 2)) - if orientation == 4: - return np.ascontiguousarray(np.flipud(arr)) - if orientation == 5: - return np.ascontiguousarray(np.swapaxes(arr, 0, 1)) - if orientation == 6: # rotate 90° CW - return np.ascontiguousarray(np.rot90(arr, 3)) - if orientation == 7: - return np.ascontiguousarray(np.rot90(np.swapaxes(arr, 0, 1), 2)) - if orientation == 8: # rotate 90° CCW - return np.ascontiguousarray(np.rot90(arr, 1)) - return arr + src = np.ascontiguousarray(arr) + h, w = src.shape[:2] + out = negpy_core.orient(src.view(np.uint8).reshape(h, w, -1), int(orientation)) + return out.view(src.dtype).reshape(out.shape[:2] + src.shape[2:]) def get_luminance(img: np.ndarray) -> np.ndarray: diff --git a/negpy_core/src/lib.rs b/negpy_core/src/lib.rs index 05508b3da..6a7b175a2 100644 --- a/negpy_core/src/lib.rs +++ b/negpy_core/src/lib.rs @@ -27,6 +27,7 @@ pub mod median; pub mod memo; pub mod meter; pub mod morph; +pub mod orient; pub mod phase; pub mod print_curve; pub mod process; diff --git a/negpy_core/src/orient.rs b/negpy_core/src/orient.rs new file mode 100644 index 000000000..13c192f78 --- /dev/null +++ b/negpy_core/src/orient.rs @@ -0,0 +1,110 @@ +//! `apply_exif_orientation`: an EXIF orientation (2-8) baked into the pixels, on any pixel size. + +use rayon::prelude::*; + +/// Output rows per parallel task, and the square tile a 90-degree turn copies at a time so +/// that both the read and the write stay in cache. +const ROWS: usize = 64; +const TILE: usize = 64; + +/// Output dimensions of an (h, w) image under `orientation`. +pub fn oriented_dims(h: usize, w: usize, orientation: u8) -> (usize, usize) { + if orientation >= 5 { (w, h) } else { (h, w) } +} + +/// The source (row, column) of output pixel (y, x), for an (h, w) source. +#[inline(always)] +fn source(orientation: u8, h: usize, w: usize, y: usize, x: usize) -> (usize, usize) { + match orientation { + 2 => (y, w - 1 - x), + 3 => (h - 1 - y, w - 1 - x), + 4 => (h - 1 - y, x), + 5 => (x, y), + 6 => (h - 1 - x, y), + 7 => (h - 1 - x, w - 1 - y), + 8 => (x, w - 1 - y), + _ => (y, x), + } +} + +fn orient_px(src: &[u8], h: usize, w: usize, orientation: u8, dst: &mut [u8]) { + let src = src.as_chunks::().0; + let (oh, ow) = oriented_dims(h, w, orientation); + dst.par_chunks_mut(ROWS * ow * N) + .enumerate() + .for_each(|(band, d)| { + let d = d.as_chunks_mut::().0; + let y0 = band * ROWS; + let rows = d.len() / ow; + for x0 in (0..ow).step_by(TILE) { + for dy in 0..rows { + let y = y0 + dy; + for x in x0..(x0 + TILE).min(ow) { + let (sy, sx) = source(orientation, h, w, y, x); + d[dy * ow + x] = src[sy * w + sx]; + } + } + } + debug_assert!(y0 + rows <= oh); + }); +} + +/// `src`, an (h, w) image of `px`-byte pixels, turned by `orientation` into `dst`. +pub fn orient(src: &[u8], h: usize, w: usize, px: usize, orientation: u8, dst: &mut [u8]) { + match px { + 1 => orient_px::<1>(src, h, w, orientation, dst), + 2 => orient_px::<2>(src, h, w, orientation, dst), + 3 => orient_px::<3>(src, h, w, orientation, dst), + 4 => orient_px::<4>(src, h, w, orientation, dst), + 6 => orient_px::<6>(src, h, w, orientation, dst), + 8 => orient_px::<8>(src, h, w, orientation, dst), + 12 => orient_px::<12>(src, h, w, orientation, dst), + 16 => orient_px::<16>(src, h, w, orientation, dst), + _ => { + let (_, ow) = oriented_dims(h, w, orientation); + dst.par_chunks_mut(ow * px).enumerate().for_each(|(y, d)| { + for (x, p) in d.chunks_exact_mut(px).enumerate() { + let (sy, sx) = source(orientation, h, w, y, x); + let o = (sy * w + sx) * px; + p.copy_from_slice(&src[o..o + px]); + } + }); + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn turns_a_2x3_image() { + // 0 1 2 + // 3 4 5 + let src: Vec = (0..6).collect(); + let run = |o, px: usize| { + let mut dst = [0u8; 6]; + orient(&src, 2, 3, px, o, &mut dst); + dst + }; + assert_eq!(run(2, 1), [2, 1, 0, 5, 4, 3]); + assert_eq!(run(3, 1), [5, 4, 3, 2, 1, 0]); + assert_eq!(run(4, 1), [3, 4, 5, 0, 1, 2]); + assert_eq!(run(5, 1), [0, 3, 1, 4, 2, 5]); + assert_eq!(run(6, 1), [3, 0, 4, 1, 5, 2]); + assert_eq!(run(7, 1), [5, 2, 4, 1, 3, 0]); + assert_eq!(run(8, 1), [2, 5, 1, 4, 0, 3]); + } + + #[test] + fn odd_pixel_sizes_take_the_generic_path() { + let src: Vec = (0..30).collect(); + let mut fast = [0u8; 18]; + let mut slow = [0u8; 30]; + orient(&src, 2, 3, 5, 6, &mut slow); + let src3: Vec = (0..18).collect(); + orient(&src3, 2, 3, 3, 6, &mut fast); + assert_eq!(&slow[..5], &src[15..20]); + assert_eq!(&fast[..3], &src3[9..12]); + } +} diff --git a/negpy_core/src/py/convert.rs b/negpy_core/src/py/convert.rs index 6a1b0b8f4..ae1a66d79 100644 --- a/negpy_core/src/py/convert.rs +++ b/negpy_core/src/py/convert.rs @@ -337,7 +337,28 @@ fn nan_to_num<'py>( }) } +/// `img`, an HxWxP uint8 view of an image with P-byte pixels, turned by EXIF `orientation`. +#[pyfunction] +fn orient<'py>( + py: Python<'py>, + img: PyReadonlyArrayDyn<'py, u8>, + orientation: u8, +) -> PyResult>> { + let [h, w, px] = img.shape()[..] else { + return Err(PyValueError::new_err("orient needs an HxWxP array")); + }; + if !(1..=8).contains(&orientation) || px == 0 { + return Err(PyValueError::new_err("orientation must be 1-8")); + } + let src = contiguous(&img)?; + let (oh, ow) = crate::orient::oriented_dims(h, w, orientation); + filled(py, &[oh, ow, px], "uint8", |d| { + crate::orient::orient(src, h, w, px, orientation, d) + }) +} + pub fn register(m: &Bound<'_, PyModule>) -> PyResult<()> { + m.add_function(wrap_pyfunction!(orient, m)?)?; m.add_function(wrap_pyfunction!(nan_to_num, m)?)?; m.add_function(wrap_pyfunction!(gray_to_rgb, m)?)?; m.add_function(wrap_pyfunction!(add, m)?)?; diff --git a/tests/test_core_orient.py b/tests/test_core_orient.py new file mode 100644 index 000000000..9c1c6ec6f --- /dev/null +++ b/tests/test_core_orient.py @@ -0,0 +1,55 @@ +import negpy_core +import numpy as np +import pytest + +from negpy.kernel.image.logic import apply_exif_orientation + +_NUMPY = { + 2: np.fliplr, + 3: lambda a: np.rot90(a, 2), + 4: np.flipud, + 5: lambda a: np.swapaxes(a, 0, 1), + 6: lambda a: np.rot90(a, 3), + 7: lambda a: np.rot90(np.swapaxes(a, 0, 1), 2), + 8: lambda a: np.rot90(a, 1), +} + + +@pytest.mark.parametrize("orientation", sorted(_NUMPY)) +@pytest.mark.parametrize( + "shape,dtype", + [ + ((67, 131, 3), np.float32), + ((131, 67), np.float32), + ((70, 90, 3), np.uint8), + ((65, 129, 3), np.uint16), + ((33, 47), np.bool_), + ((9, 11, 5), np.uint8), + ], +) +def test_orientation_matches_numpy(orientation, shape, dtype): + rng = np.random.default_rng(orientation) + arr = (rng.random(shape) * 200).astype(dtype) + out = apply_exif_orientation(arr, orientation) + expected = _NUMPY[orientation](arr) + assert out.dtype == arr.dtype and out.shape == expected.shape + assert out.flags.c_contiguous + assert np.array_equal(out, expected) + + +def test_orientation_of_a_strided_view(): + arr = np.arange(40 * 60 * 3, dtype=np.float32).reshape(40, 60, 3)[5:30, 10:50] + assert np.array_equal(apply_exif_orientation(arr, 6), np.rot90(arr, 3)) + + +def test_orient_rejects_bad_input(): + with pytest.raises(ValueError): + negpy_core.orient(np.zeros((2, 2, 1), np.uint8), 9) + with pytest.raises(ValueError): + negpy_core.orient(np.zeros((2, 2), np.uint8), 6) + + +def test_identity_and_unknown_orientations_return_the_input(): + arr = np.zeros((3, 4), np.float32) + for orientation in (None, 0, 1, 9): + assert apply_exif_orientation(arr, orientation) is arr From 38bfcfebcff6f08f46eee71f3a51430838363ba3 Mon Sep 17 00:00:00 2001 From: Marcin Zawalski Date: Sun, 25 Oct 2026 09:11:46 +0100 Subject: [PATCH 02/12] perf(core): IR downsample rows and IR scaling in Rust _erode_resize_bounded ran its vertical INTER_AREA pass as a Python loop of small NumPy calls per output row. negpy_core.area_rows_into runs the same float64 weights and float32 products, summed down the rows in NumPy's order, in parallel; the output is bit-identical (NumPy oracle test). normalize_ir_to_float32 scales uint8/uint16 IR planes in one pass. 23.5 MP IR plane, 24 threads: downsample_ir to 2000 px 68 -> 23 ms, normalize 11.4 -> 9.5 ms. Both run on every IR frame open. --- .github/workflows/ci.yml | 2 +- negpy/features/retouch/logic.py | 13 +---- negpy/infrastructure/loaders/ir_planes.py | 7 ++- negpy_core/src/convert.rs | 6 +++ negpy_core/src/py/convert.rs | 24 +++++++++ negpy_core/src/py/resize.rs | 39 ++++++++++++++ negpy_core/src/resize.rs | 41 +++++++++++++++ tests/test_core_ir_downsample.py | 62 +++++++++++++++++++++++ 8 files changed, 178 insertions(+), 16 deletions(-) create mode 100644 tests/test_core_ir_downsample.py diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 33f7cdc5a..799079d61 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -129,4 +129,4 @@ jobs: run: uv sync --group dev --group plustek - name: negpy_core pins and oracles - run: uv run pytest tests/test_percentiles.py tests/test_core_same_bits.py tests/test_core_orient.py tests/test_core_resize.py tests/test_core_filter.py tests/test_core_morph.py tests/test_core_gray.py tests/test_core_contours.py tests/test_core_shapes.py tests/test_core_components.py tests/test_core_color.py tests/test_core_warp.py tests/test_core_lab.py tests/test_core_altprocess.py + run: uv run pytest tests/test_percentiles.py tests/test_core_same_bits.py tests/test_core_orient.py tests/test_core_ir_downsample.py tests/test_core_resize.py tests/test_core_filter.py tests/test_core_morph.py tests/test_core_gray.py tests/test_core_contours.py tests/test_core_shapes.py tests/test_core_components.py tests/test_core_color.py tests/test_core_warp.py tests/test_core_lab.py tests/test_core_altprocess.py diff --git a/negpy/features/retouch/logic.py b/negpy/features/retouch/logic.py index b50657829..cd185df9b 100644 --- a/negpy/features/retouch/logic.py +++ b/negpy/features/retouch/logic.py @@ -674,17 +674,8 @@ def _erode_resize_bounded(plane: np.ndarray, dims: Tuple[int, int], kernel: np.n eroded = negpy_core.morphology_ex(source, negpy_core.MORPH_ERODE, kernel) core = np.ascontiguousarray(eroded[source_top - read_top : source_bottom - read_top]) horizontal = negpy_core.resize_area(core, (dw, len(core))) - # Keep the full-image sampling grid across fractional block boundaries. - for row in range(top, bottom): - start = row * scale_y - end = start + scale_y - first, last = math.floor(start), min(h, math.ceil(end)) - indices = np.arange(first, last) - overlap = np.minimum(indices + 1, end) - np.maximum(indices, start) - # OpenCV's area table omits fractional edges at or below this cutoff. - weights = (np.where(overlap > 1e-3, overlap, 0.0) / min(scale_y, h - start)).astype(np.float32) - values = horizontal[first - source_top : last - source_top] - output[row] = np.sum(values * weights.reshape((-1,) + (1,) * (values.ndim - 1)), axis=0) + # The vertical pass keeps the full-image sampling grid across fractional block boundaries. + negpy_core.area_rows_into(horizontal, output[top:bottom], top, source_top, scale_y, h) return output diff --git a/negpy/infrastructure/loaders/ir_planes.py b/negpy/infrastructure/loaders/ir_planes.py index 9cdeb648f..c5a7c5b8b 100644 --- a/negpy/infrastructure/loaders/ir_planes.py +++ b/negpy/infrastructure/loaders/ir_planes.py @@ -8,6 +8,7 @@ from typing import Any, Iterable, Optional +import negpy_core import numpy as np _NEW_SUBFILE_TYPE = 254 @@ -19,10 +20,8 @@ def normalize_ir_to_float32(ir: np.ndarray) -> np.ndarray: """Single-channel uint8/uint16/float ndarray → float32 in [0,1].""" if ir.ndim == 3: ir = ir[:, :, 0] - if ir.dtype == np.uint8: - return ir.astype(np.float32) * (1.0 / 255.0) - if ir.dtype == np.uint16: - return ir.astype(np.float32) * (1.0 / 65535.0) + if ir.dtype in (np.uint8, np.uint16): + return negpy_core.int_to_unit_float32(np.ascontiguousarray(ir)) return np.clip(ir.astype(np.float32), 0.0, 1.0) diff --git a/negpy_core/src/convert.rs b/negpy_core/src/convert.rs index 780043b43..120f0ad6b 100644 --- a/negpy_core/src/convert.rs +++ b/negpy_core/src/convert.rs @@ -41,6 +41,12 @@ pub fn uint16_to_float32(src: &[u16], dst: &mut [f32]) { crate::par_zip(dst, src, move |d, &v| *d = (v as f32 as f64 * inv) as f32); } +/// `v.astype(np.float32) * scale` with `scale` already float32, NumPy's NEP 50 rule for a +/// Python float times a float32 array. +pub fn scaled_to_float32 + Sync>(src: &[T], scale: f32, dst: &mut [f32]) { + crate::par_zip(dst, src, move |d, &v| *d = v.into() * scale); +} + /// Numba's fastmath luminance in float64: two fused multiply-adds. fn luma(p: &[f32]) -> f64 { LUMA_B.mul_add( diff --git a/negpy_core/src/py/convert.rs b/negpy_core/src/py/convert.rs index ae1a66d79..58b33a68c 100644 --- a/negpy_core/src/py/convert.rs +++ b/negpy_core/src/py/convert.rs @@ -337,6 +337,29 @@ fn nan_to_num<'py>( }) } +/// `ir.astype(np.float32) * (1 / 255)` or `* (1 / 65535)` of a contiguous uint8 or uint16 array. +#[pyfunction] +fn int_to_unit_float32<'py>( + py: Python<'py>, + img: &Bound<'py, PyAny>, +) -> PyResult> { + if let Ok(a) = img.extract::>() { + let src = contiguous(&a)?; + let scale = (1.0f64 / 65535.0) as f32; + return Ok(filled(py, a.shape(), "float32", |d| { + convert::scaled_to_float32(src, scale, d) + })? + .into_any()); + } + let a = img.extract::>()?; + let src = contiguous(&a)?; + let scale = (1.0f64 / 255.0) as f32; + Ok(filled(py, a.shape(), "float32", |d| { + convert::scaled_to_float32(src, scale, d) + })? + .into_any()) +} + /// `img`, an HxWxP uint8 view of an image with P-byte pixels, turned by EXIF `orientation`. #[pyfunction] fn orient<'py>( @@ -359,6 +382,7 @@ fn orient<'py>( pub fn register(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_function(wrap_pyfunction!(orient, m)?)?; + m.add_function(wrap_pyfunction!(int_to_unit_float32, m)?)?; m.add_function(wrap_pyfunction!(nan_to_num, m)?)?; m.add_function(wrap_pyfunction!(gray_to_rgb, m)?)?; m.add_function(wrap_pyfunction!(add, m)?)?; diff --git a/negpy_core/src/py/resize.rs b/negpy_core/src/py/resize.rs index 2d4f69aea..6eb46078b 100644 --- a/negpy_core/src/py/resize.rs +++ b/negpy_core/src/py/resize.rs @@ -298,7 +298,46 @@ fn resize_rows_linear<'py>( } } +/// `area_rows_into(src, dst, top, source_top, scale_y, full_h)`: the vertical INTER_AREA rows +/// `top..top + len(dst)` of `_erode_resize_bounded`, from a block's horizontally resized rows. +#[pyfunction] +fn area_rows_into( + py: Python<'_>, + src: PyReadonlyArrayDyn<'_, f32>, + mut dst: numpy::PyReadwriteArrayDyn<'_, f32>, + top: usize, + source_top: usize, + scale_y: f64, + full_h: usize, +) -> PyResult<()> { + let len: usize = src.shape()[1..].iter().product(); + if src.ndim() < 2 || dst.shape()[1..].iter().product::() != len || len == 0 { + return Err(PyValueError::new_err( + "needs src and dst rows of one length", + )); + } + let rows = src.shape()[0]; + let n = dst.shape()[0]; + let first_needed = (top as f64 * scale_y).floor() as usize; + let last_needed = + (((top + n).saturating_sub(1) as f64 * scale_y + scale_y).ceil() as usize).min(full_h); + if n > 0 && (first_needed < source_top || last_needed > source_top + rows) { + return Err(PyValueError::new_err("rows outside the block")); + } + let s = src + .as_slice() + .map_err(|_| PyValueError::new_err("needs a contiguous src"))?; + let d = dst + .as_slice_mut() + .map_err(|_| PyValueError::new_err("needs a contiguous dst"))?; + super::detach(py, d.len(), || { + resize::area_rows(s, len, source_top, top, scale_y, full_h, d) + }); + Ok(()) +} + pub fn register(m: &Bound<'_, PyModule>) -> PyResult<()> { + m.add_function(wrap_pyfunction!(area_rows_into, m)?)?; m.add_function(wrap_pyfunction!(resize_rows_linear, m)?)?; m.add_function(wrap_pyfunction!(resize_u8, m)?)?; m.add_function(wrap_pyfunction!(resize_kernel, m)?)?; diff --git a/negpy_core/src/resize.rs b/negpy_core/src/resize.rs index 5ddceed43..ea8a4ffac 100644 --- a/negpy_core/src/resize.rs +++ b/negpy_core/src/resize.rs @@ -1010,6 +1010,47 @@ pub fn expand_plane( Some(()) } +/// The vertical INTER_AREA pass of `retouch.logic._erode_resize_bounded` on one block: output +/// rows `top..top + dst rows` of a `full_h`-row image shrunk by `scale_y`, from the block's +/// horizontally resized rows `src` (row `source_top` first, `len` floats per row). Weights in +/// float64 then float32; each product in float32, summed down the rows in order, as NumPy does. +pub fn area_rows( + src: &[f32], + len: usize, + source_top: usize, + top: usize, + scale_y: f64, + full_h: usize, + dst: &mut [f32], +) { + use rayon::prelude::*; + let src_rows = src.len() / len; + dst.par_chunks_mut(len).enumerate().for_each(|(r, out)| { + let start = (top + r) as f64 * scale_y; + let end = start + scale_y; + let first = start.floor() as usize; + let last = (end.ceil() as usize).min(full_h); + let norm = scale_y.min(full_h as f64 - start); + out.fill(0.0); + for (k, i) in (first..last).enumerate() { + let overlap = ((i + 1) as f64).min(end) - (i as f64).max(start); + let w = ((if overlap > 1e-3 { overlap } else { 0.0 }) / norm) as f32; + let row = i - source_top; + debug_assert!(row < src_rows); + let vals = &src[row * len..(row + 1) * len]; + if k == 0 { + for (o, &v) in out.iter_mut().zip(vals) { + *o = v * w; + } + } else { + for (o, &v) in out.iter_mut().zip(vals) { + *o += v * w; + } + } + } + }); +} + #[cfg(test)] mod tests { use super::*; diff --git a/tests/test_core_ir_downsample.py b/tests/test_core_ir_downsample.py new file mode 100644 index 000000000..d299728fc --- /dev/null +++ b/tests/test_core_ir_downsample.py @@ -0,0 +1,62 @@ +import math + +import negpy_core +import numpy as np +import pytest + +import negpy.features.retouch.logic as retouch +from negpy.infrastructure.loaders.ir_planes import normalize_ir_to_float32 + + +def _numpy_rows(horizontal, output, top, bottom, source_top, scale_y, h): + """The NumPy vertical pass negpy_core.area_rows_into replaces, as the oracle.""" + for row in range(top, bottom): + start = row * scale_y + end = start + scale_y + first, last = math.floor(start), min(h, math.ceil(end)) + indices = np.arange(first, last) + overlap = np.minimum(indices + 1, end) - np.maximum(indices, start) + weights = (np.where(overlap > 1e-3, overlap, 0.0) / min(scale_y, h - start)).astype(np.float32) + values = horizontal[first - source_top : last - source_top] + output[row] = np.sum(values * weights.reshape((-1,) + (1,) * (values.ndim - 1)), axis=0) + + +@pytest.mark.parametrize("shape,dims", [((997, 641), (233, 150)), ((1200, 800, 3), (301, 200)), ((640, 999), (128, 81))]) +def test_area_rows_match_numpy(shape, dims): + rng = np.random.default_rng(len(shape) + dims[0]) + plane = rng.random(shape, dtype=np.float32) + h = shape[0] + dw, dh = dims + scale_y = h / dh + for top, bottom in [(0, dh // 3), (dh // 3, dh - 7), (dh - 7, dh)]: + source_top = math.floor(top * scale_y) + source_bottom = min(h, math.ceil((bottom - 1) * scale_y + scale_y)) + horizontal = negpy_core.resize_area(np.ascontiguousarray(plane[source_top:source_bottom]), (dw, source_bottom - source_top)) + expected = np.zeros((dh, dw) + shape[2:], np.float32) + got = np.zeros_like(expected) + _numpy_rows(horizontal, expected, top, bottom, source_top, scale_y, h) + negpy_core.area_rows_into(horizontal, got[top:bottom], top, source_top, scale_y, h) + assert np.array_equal(got[top:bottom], expected[top:bottom]) + + +def test_bounded_downsample_matches_whole_plane(monkeypatch): + rng = np.random.default_rng(7) + plane = rng.random((1500, 1001), dtype=np.float32) + whole = retouch.downsample_ir(plane, 400) + monkeypatch.setattr(retouch, "_IR_DOWNSAMPLE_WORK_BYTES", 1 << 18) + blocks = retouch.downsample_ir(plane, 400) + assert blocks.shape == whole.shape + assert np.allclose(blocks, whole, atol=1e-6) + + +def test_area_rows_rejects_rows_outside_the_block(): + src = np.zeros((4, 10), np.float32) + with pytest.raises(ValueError): + negpy_core.area_rows_into(src, np.zeros((5, 10), np.float32), 0, 0, 3.0, 30) + + +@pytest.mark.parametrize("dtype,scale", [(np.uint16, 1.0 / 65535.0), (np.uint8, 1.0 / 255.0)]) +def test_ir_normalize_matches_numpy(dtype, scale): + ir = np.random.default_rng(3).integers(0, np.iinfo(dtype).max + 1, (311, 257), dtype=dtype) + assert np.array_equal(normalize_ir_to_float32(ir), ir.astype(np.float32) * scale) + assert np.array_equal(normalize_ir_to_float32(ir[:, :, None].repeat(2, axis=2)), ir.astype(np.float32) * scale) From b53233869f238999d865ce4082b613f2c5e83b0e Mon Sep 17 00:00:00 2001 From: Marcin Zawalski Date: Sun, 25 Oct 2026 09:13:57 +0100 Subject: [PATCH 03/12] perf(export): 8-bit export without PIL round trips _export_pixels sent the quantized buffer through Image.fromarray and np.asarray twice around the Rust matrix/TRC kernel, and the kernel's decode/encode tables were rebuilt from the ICC bytes on every export. apply_color_management now takes and returns a uint8 array (PIL only for the lcms fallback) and the tables are cached per profile pair. Pixels and embedded profile are unchanged. 23.5 MP, Adobe RGB -> sRGB: _export_pixels 173 -> 49 ms per JPEG/WebP export. --- negpy/infrastructure/display/icc_lut.py | 17 ++++++++-- negpy/services/rendering/image_processor.py | 36 +++++++++------------ tests/test_image_processor.py | 6 ++-- 3 files changed, 33 insertions(+), 26 deletions(-) diff --git a/negpy/infrastructure/display/icc_lut.py b/negpy/infrastructure/display/icc_lut.py index 7cb1ce18e..7864d6a60 100644 --- a/negpy/infrastructure/display/icc_lut.py +++ b/negpy/infrastructure/display/icc_lut.py @@ -8,7 +8,8 @@ per-sample fidelity is 8-bit (bounded by LUT size and profile non-linearity). """ -from typing import Any, Optional +from functools import lru_cache +from typing import Any, Optional, Tuple import negpy_core import numpy as np @@ -101,6 +102,16 @@ def apply_matrix_trc_u8(img_u8: np.ndarray, src_icc: bytes, dst_icc: bytes) -> O """Relative-colorimetric transform between two matrix/TRC profiles on a (H,W,3) uint8 image: TRC decode, inv(dst_pcs) @ src_pcs with gamut clip, TRC encode. None when either profile is not matrix/TRC; the caller falls back to lcms.""" + tables = _matrix_trc_tables(src_icc, dst_icc) + if tables is None: + return None + dec_lut, m32, enc_lut = tables + return negpy_core.matrix_trc_u8(np.ascontiguousarray(img_u8), dec_lut, m32, enc_lut) + + +@lru_cache(maxsize=16) +def _matrix_trc_tables(src_icc: bytes, dst_icc: bytes) -> Optional[Tuple[np.ndarray, np.ndarray, np.ndarray]]: + """(decode LUT, matrix, encode LUT) for a matrix/TRC profile pair, or None.""" from negpy.infrastructure.display.icc_profile import ( extract_pcs_matrix, extract_trc_decode_samples, @@ -131,4 +142,6 @@ def apply_matrix_trc_u8(img_u8: np.ndarray, src_icc: bytes, dst_icc: bytes) -> O enc_lut = np.ascontiguousarray(np.stack(enc_rows).astype(np.float32)) dec_lut = np.ascontiguousarray(np.clip(dec, 0.0, 1.0).astype(np.float32)) m32 = np.ascontiguousarray(m.astype(np.float32)) - return negpy_core.matrix_trc_u8(np.ascontiguousarray(img_u8), dec_lut, m32, enc_lut) + for table in (dec_lut, m32, enc_lut): + table.setflags(write=False) + return dec_lut, m32, enc_lut diff --git a/negpy/services/rendering/image_processor.py b/negpy/services/rendering/image_processor.py index 10de509a1..177a1a5c9 100644 --- a/negpy/services/rendering/image_processor.py +++ b/negpy/services/rendering/image_processor.py @@ -1853,11 +1853,8 @@ def _export_pixels( return self._apply_color_management_u16(float_to_uint16(buffer), working_color_space, color_space, icc_output, icc_input) if is_greyscale: pil_img, icc_bytes = self._greyscale_to_pil_u8(buffer, working_color_space, color_space, icc_output, icc_input) - else: - pil_img, icc_bytes = self.apply_color_management( - Image.fromarray(float_to_uint8(buffer)), working_color_space, color_space, icc_output, icc_input - ) - return np.asarray(pil_img), icc_bytes + return np.asarray(pil_img), icc_bytes + return self.apply_color_management(float_to_uint8(buffer), working_color_space, color_space, icc_output, icc_input) def render_display_array( self, @@ -2231,37 +2228,35 @@ def _log_cms_codec_dlopen_error() -> None: def apply_color_management( self, - pil_img: Image.Image, + img: np.ndarray, working_color_space: str, color_space: str, output_icc_path: Optional[str], input_icc_path: Optional[str] = None, - ) -> Tuple[Image.Image, Optional[bytes]]: - """ICC transform for export. Source is the input override or working space; - destination is the output override or target space.""" + ) -> Tuple[np.ndarray, Optional[bytes]]: + """ICC transform of a uint8 HxWx3 or HxW export buffer. Source is the input override or + working space; destination is the output override or target space.""" has_custom = self._has_custom_icc(input_icc_path, output_icc_path) if not has_custom and working_color_space == color_space: - return pil_img, self._get_target_icc_bytes(color_space, None) + return img, self._get_target_icc_bytes(color_space, None) try: p_src = self._resolve_src_profile(working_color_space, input_icc_path) p_dst = self._resolve_dst_profile(color_space, output_icc_path) if p_dst is None: - return pil_img, self._get_target_icc_bytes(color_space, None) - - if pil_img.mode not in ("RGB", "L"): - pil_img = pil_img.convert("RGB" if pil_img.mode != "I;16" else "L") + return img, self._get_target_icc_bytes(color_space, None) - if pil_img.mode == "RGB": + if img.ndim == 3: # Matrix/TRC pairs take the parallel kernel; LUT profiles (printer # ICCs) return None and fall through to exact lcms. src_bytes = self._get_target_icc_bytes(working_color_space, input_icc_path) dst_bytes = self._get_target_icc_bytes(color_space, output_icc_path) if src_bytes and dst_bytes: - fast = apply_matrix_trc_u8(np.asarray(pil_img), src_bytes, dst_bytes) + fast = apply_matrix_trc_u8(img, src_bytes, dst_bytes) if fast is not None: - return Image.fromarray(fast), dst_bytes + return fast, dst_bytes + pil_img = Image.fromarray(img) result_pil = ImageCms.profileToProfile( pil_img, p_src, @@ -2270,14 +2265,13 @@ def apply_color_management( outputMode="RGB" if pil_img.mode != "L" else "L", flags=ImageCms.Flags.BLACKPOINTCOMPENSATION, ) - if result_pil: - pil_img = result_pil - return pil_img, self._get_target_icc_bytes(color_space, output_icc_path) + out = np.asarray(result_pil) if result_pil else img + return out, self._get_target_icc_bytes(color_space, output_icc_path) except Exception as e: logger.error(f"CMS transformation failed: {e}") # The pixels never left the working space, so tag them with it (a CMYK or Gray # profile cannot take RGB pixels). Untagged, the export reads as sRGB. - return pil_img, self._get_target_icc_bytes(working_color_space, input_icc_path) + return img, self._get_target_icc_bytes(working_color_space, input_icc_path) @staticmethod def soft_proof_preview( diff --git a/tests/test_image_processor.py b/tests/test_image_processor.py index 99ce51cea..f1bd36e39 100644 --- a/tests/test_image_processor.py +++ b/tests/test_image_processor.py @@ -463,16 +463,16 @@ def test_detect_luma_reads_the_prepared_plane(monkeypatch) -> None: def test_an_8bit_export_to_a_profile_lcms_refuses_stays_tagged_with_the_working_space(tmp_path): """A CMYK, Gray or Lab Export ICC cannot take RGB pixels; untagged, they read as sRGB.""" - from PIL import Image, ImageCms + from PIL import ImageCms from negpy.infrastructure.display.color_spaces import WORKING_COLOR_SPACE path = tmp_path / "lab.icc" path.write_bytes(ImageCms.ImageCmsProfile(ImageCms.createProfile("LAB")).tobytes()) - img = Image.fromarray((np.random.default_rng(0).random((8, 8, 3)) * 255).astype(np.uint8)) + img = (np.random.default_rng(0).random((8, 8, 3)) * 255).astype(np.uint8) proc = ImageProcessor() out, icc = proc.apply_color_management(img, WORKING_COLOR_SPACE, WORKING_COLOR_SPACE, str(path)) - assert out.mode == "RGB" + assert out.shape == (8, 8, 3) and out.dtype == np.uint8 assert icc == proc._get_target_icc_bytes(WORKING_COLOR_SPACE, None) From 857aab3089b109b398ac0ba504a73264a513c825 Mon Sep 17 00:00:00 2001 From: Marcin Zawalski Date: Sun, 25 Oct 2026 09:18:51 +0100 Subject: [PATCH 04/12] fix(export): exact 16-bit color transform in Rust 16-bit RGB exports (TIFF, PNG, JXL) went through lcms, which resamples a 16-bit transform into a grid table. Against lcms with cmsFLAGS_NOOPTIMIZE that grid is off by up to 5386 levels (8%) for Adobe RGB -> sRGB and 3764 for P3 near the destination gamut edge, mean 25. 8-bit exports already use the exact matrix/TRC kernel, so 16-bit exports were the less accurate ones. matrix_trc_u16 mirrors matrix_trc_u8: per-code TRC decode, the PCS matrix with a gamut clip, and an encode table sampled at linear t^2 so the steep end near black interpolates finely. It is within one level of exact lcms on every bundled export space. LUT profiles keep lcms. 23.5 MP, Adobe RGB -> sRGB: 124 -> 43 ms per 16-bit export. 16-bit export pixels change by the lcms grid error above; thumbnails are not affected. --- docs/PIPELINE.md | 1 + negpy/infrastructure/display/icc_lut.py | 51 ++++++++++++++ negpy/services/rendering/image_processor.py | 9 ++- negpy_core/src/display.rs | 73 +++++++++++++++++++++ negpy_core/src/py/display.rs | 34 ++++++++++ tests/test_core_same_bits.py | 4 ++ tests/test_export_color_management.py | 34 ++++++++-- 7 files changed, 199 insertions(+), 7 deletions(-) diff --git a/docs/PIPELINE.md b/docs/PIPELINE.md index e766e4281..8f51603cd 100644 --- a/docs/PIPELINE.md +++ b/docs/PIPELINE.md @@ -99,6 +99,7 @@ Pixel values are linear transmittance. The source corrections run in that domain * **Waveform**: encoded luma in 128 bins for each of 256 columns (`WAVEFORM_BINS`, `WAVEFORM_COLS`), with each cell's mean color (`WAVEFORM_COLOR_SCALE`) and a per-channel parade. The density pass bins scan density per crop column for the scan Parade (`SCAN_DENSITY_RANGE`, inverting the stretch) and counts linear pixels at `scan_clip_level` before the unmix. The base line is `normalization.base_scan_density`. * **Vectorscope**: CIELAB a\*b\* in 180 hue × 64 chroma bins, log chroma with a knee at 4, last bin past chroma 60 (`SCOPE_*`). * **Input ICC**: overrides primaries only. A matrix/TRC profile contributes its primaries (Bradford-adapted to D65 through its `chad` tag) as a plain matrix, so its TRC is inert. A LUT (A2B0/B2A0) profile runs the full CMS transform. + * **Export transform**: between two matrix/TRC profiles (every bundled export space), 8- and 16-bit exports run `matrix_trc_u8`/`matrix_trc_u16`: TRC decode, $M_{dst,PCS}^{-1} M_{src,PCS}$ with a gamut clip, TRC encode. The 16-bit encode table samples linear values at $t^2$, so it stays within one level of lcms without lcms's grid optimization. LUT profiles go through lcms. ### Automatic helpers diff --git a/negpy/infrastructure/display/icc_lut.py b/negpy/infrastructure/display/icc_lut.py index 7864d6a60..dd4f15392 100644 --- a/negpy/infrastructure/display/icc_lut.py +++ b/negpy/infrastructure/display/icc_lut.py @@ -145,3 +145,54 @@ def _matrix_trc_tables(src_icc: bytes, dst_icc: bytes) -> Optional[Tuple[np.ndar for table in (dec_lut, m32, enc_lut): table.setflags(write=False) return dec_lut, m32, enc_lut + + +_ENC16_SIZE = 16385 +_DEC_DENSE = 1 << 20 + + +def apply_matrix_trc_u16(img_u16: np.ndarray, src_icc: bytes, dst_icc: bytes) -> Optional[np.ndarray]: + """The 16-bit twin of ``apply_matrix_trc_u8``: exact per-code TRC decode, the PCS matrix + with a gamut clip, and an encode table fine enough for 16 bits. None when either profile + is not matrix/TRC; the caller falls back to lcms.""" + tables = _matrix_trc_tables_u16(src_icc, dst_icc) + if tables is None: + return None + dec_lut, m32, enc_lut = tables + return negpy_core.matrix_trc_u16(np.ascontiguousarray(img_u16), dec_lut, m32, enc_lut) + + +@lru_cache(maxsize=16) +def _matrix_trc_tables_u16(src_icc: bytes, dst_icc: bytes) -> Optional[Tuple[np.ndarray, np.ndarray, np.ndarray]]: + """(3 x 65536 decode, matrix, 3 x K encode at linear t^2) for a matrix/TRC pair, or None.""" + from negpy.infrastructure.display.icc_profile import ( + extract_pcs_matrix, + extract_trc_decode_samples, + is_matrix_trc_profile, + ) + + if not (is_matrix_trc_profile(src_icc) and is_matrix_trc_profile(dst_icc)): + return None + m_src = extract_pcs_matrix(src_icc) + m_dst = extract_pcs_matrix(dst_icc) + if m_src is None or m_dst is None: + return None + dec = extract_trc_decode_samples(src_icc, np.arange(65536) / 65535.0) + x_enc = np.linspace(0.0, 1.0, _DEC_DENSE) + dst_dec = extract_trc_decode_samples(dst_icc, x_enc) + if dec is None or dst_dec is None: + return None + try: + m = np.linalg.inv(m_dst) @ m_src + except np.linalg.LinAlgError: + return None + t = np.linspace(0.0, 1.0, _ENC16_SIZE) + enc_rows = [np.interp(t * t, np.maximum.accumulate(np.clip(dst_dec[ch], 0.0, 1.0)), x_enc) for ch in range(3)] + tables = ( + np.ascontiguousarray(np.clip(dec, 0.0, 1.0).astype(np.float32)), + np.ascontiguousarray(m.astype(np.float32)), + np.ascontiguousarray(np.stack(enc_rows).astype(np.float32)), + ) + for table in tables: + table.setflags(write=False) + return tables diff --git a/negpy/services/rendering/image_processor.py b/negpy/services/rendering/image_processor.py index 177a1a5c9..9a05f66cc 100644 --- a/negpy/services/rendering/image_processor.py +++ b/negpy/services/rendering/image_processor.py @@ -104,7 +104,7 @@ from negpy.services.export.print import PrintService from negpy.services.export.encoders import encode_jpeg, encode_jxl, encode_png, encode_tiff, encode_webp from negpy.infrastructure.display.color_spaces import ColorSpaceRegistry, WORKING_COLOR_SPACE -from negpy.infrastructure.display.icc_lut import DEFAULT_LUT_SIZE, apply_icc_u16_rgb, apply_matrix_trc_u8 +from negpy.infrastructure.display.icc_lut import DEFAULT_LUT_SIZE, apply_icc_u16_rgb, apply_matrix_trc_u16, apply_matrix_trc_u8 # Preview soft-proof LUT grid. Finer than the display LUT because the proof clips at the # output gamut boundary, and interpolating across that kink is where the error is. @@ -2188,6 +2188,13 @@ def _apply_color_management_u16( logger.warning("CMS skipped: ICC profile not found") return img_u16, self._get_target_icc_bytes(color_space, output_icc_path) + # Matrix/TRC pairs take the exact kernel: lcms resamples a 16-bit transform to a + # grid, off by up to 8% near the destination gamut edge. LUT profiles keep lcms. + if img_u16.ndim == 3: + fast = apply_matrix_trc_u16(img_u16, src_bytes, dst_bytes) + if fast is not None: + return fast, dst_bytes + try: result = _cms_transform_strips(np.ascontiguousarray(img_u16), src_bytes, dst_bytes) except ImportError as e: diff --git a/negpy_core/src/display.rs b/negpy_core/src/display.rs index e86ee3fe8..588fdc3dd 100644 --- a/negpy_core/src/display.rs +++ b/negpy_core/src/display.rs @@ -171,3 +171,76 @@ fn matrix_trc_block(img: &[u8], out: &mut [u8], dec: &[f32], m: [[f32; 3]; 3], e } } } + +/// The 16-bit matrix/TRC transform: TRC decode through `dec` (3 x 65536, one entry per code), +/// the 3x3 `m` with a gamut clip, TRC encode through `enc` (3 x K) sampled at linear values +/// t^2 for t evenly spaced in 0..1, so a steep encode near black still interpolates finely; +/// to 16-bit rounded half up. +pub fn matrix_trc_u16(img: &[u16], dec: &[f32], m: [[f32; 3]; 3], enc: &[f32], out: &mut [u16]) { + out.par_chunks_mut(3 * 4096) + .zip(img.par_chunks(3 * 4096)) + .for_each(|(o, i)| matrix_trc_u16_block(i, o, dec, m, enc)); +} + +#[multiversion(targets( + "x86_64+avx512f+avx512bw+avx512vl+avx512dq", + "x86_64+avx2+fma", + "aarch64+neon" +))] +fn matrix_trc_u16_block(img: &[u16], out: &mut [u16], dec: &[f32], m: [[f32; 3]; 3], enc: &[f32]) { + let k = enc.len() / 3; + let kmax = (k - 1) as f64; + let mut rgb = [[0f32; 3]; RUN]; + let (mut idx, mut frac) = ([0usize; RUN], [0f64; RUN]); + for (o, p) in out.chunks_mut(3 * RUN).zip(img.chunks(3 * RUN)) { + let (o, p) = (o.as_chunks_mut::<3>().0, p.as_chunks::<3>().0); + let n = p.len(); + for (q, p) in rgb.iter_mut().zip(p) { + *q = std::array::from_fn(|c| dec[c * 65536 + p[c] as usize]); + } + for c in 0..3 { + let row = m[c]; + for i in 0..n { + let q = rgb[i]; + let v = q[2] + .mul_add(row[2], q[1].mul_add(row[1], row[0] * q[0])) + .clamp(0.0, 1.0); + let f = (v as f64).sqrt() * kmax; + let i0 = (crate::convert::trunc_bits(f) as usize).min(k - 2); + idx[i] = c * k + i0; + frac[i] = f - i0 as f64; + } + for i in 0..n { + let e = lerp(enc[idx[i]] as f64, enc[idx[i] + 1] as f64, frac[i]); + o[i][c] = + crate::convert::trunc_bits((e * 65535.0 + 0.5).clamp(0.0, 65535.0)) as u16; + } + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn identity_u16_round_trips_every_code() { + let k = 4097; + let dec: Vec = (0..3 * 65536) + .map(|i| (i % 65536) as f32 / 65535.0) + .collect(); + let enc: Vec = (0..3 * k) + .map(|i| { + let t = (i % k) as f64 / (k - 1) as f64; + (t * t) as f32 + }) + .collect(); + let eye = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]; + let img: Vec = (0..3 * 30000).map(|i| ((i * 7) % 65536) as u16).collect(); + let mut out = vec![0u16; img.len()]; + matrix_trc_u16(&img, &dec, eye, &enc, &mut out); + for (a, b) in img.iter().zip(&out) { + assert!((*a as i32 - *b as i32).abs() <= 1, "{a} -> {b}"); + } + } +} diff --git a/negpy_core/src/py/display.rs b/negpy_core/src/py/display.rs index 693723ace..ee076d36f 100644 --- a/negpy_core/src/py/display.rs +++ b/negpy_core/src/py/display.rs @@ -99,7 +99,41 @@ fn matrix_trc_u8<'py>( }) } +/// `matrix_trc_u16(img, dec_lut, m, enc_lut)`: the 16-bit matrix/TRC transform of an HxWx3 +/// uint16 image, a 3 x 65536 decode, a 3x3 matrix and a 3 x K encode sampled at t^2. +#[pyfunction] +fn matrix_trc_u16<'py>( + py: Python<'py>, + img: PyReadonlyArrayDyn<'py, u16>, + dec_lut: PyReadonlyArray2<'py, f32>, + m: PyReadonlyArray2<'py, f32>, + enc_lut: PyReadonlyArray2<'py, f32>, +) -> PyResult>> { + let src = rgb(&img)?; + if dec_lut.shape() != [3, 65536] + || m.shape() != [3, 3] + || enc_lut.shape()[0] != 3 + || enc_lut.shape()[1] < 2 + { + return Err(PyValueError::new_err( + "needs a 3 x 65536 decode, a 3x3 matrix and a 3 x K encode", + )); + } + let dec = dec_lut + .as_slice() + .map_err(|_| PyValueError::new_err("needs contiguous tables"))?; + let enc = enc_lut + .as_slice() + .map_err(|_| PyValueError::new_err("needs contiguous tables"))?; + let m = m.as_array(); + let m = std::array::from_fn(|r| std::array::from_fn(|c| m[[r, c]])); + filled(py, img.shape(), |d| { + display::matrix_trc_u16(src, dec, m, enc, d) + }) +} + pub fn register(m: &Bound<'_, PyModule>) -> PyResult<()> { + m.add_function(wrap_pyfunction!(matrix_trc_u16, m)?)?; m.add_function(wrap_pyfunction!(apply_lut_u16, m)?)?; m.add_function(wrap_pyfunction!(apply_lut_f32, m)?)?; m.add_function(wrap_pyfunction!(matrix_trc_u8, m)?) diff --git a/tests/test_core_same_bits.py b/tests/test_core_same_bits.py index c5accec89..3a94471d0 100644 --- a/tests/test_core_same_bits.py +++ b/tests/test_core_same_bits.py @@ -73,6 +73,9 @@ def _filled(fill): ), "lut_f32": lambda: negpy_core.apply_lut_f32(_IMG3, _LUT), "matrix_trc_u8": lambda: negpy_core.matrix_trc_u8((_IMG3 * 255).astype(np.uint8), np.ascontiguousarray(_RAMP[:, ::256]), _M, _RAMP), + "matrix_trc_u16": lambda: negpy_core.matrix_trc_u16( + (_IMG3 * 65535).astype(np.uint16), np.ascontiguousarray(np.tile(_LIN, (3, 1))), _M, np.ascontiguousarray(_RAMP[:, ::16]) + ), "usm_sharpen": lambda: negpy_core.usm_sharpen( _IMG3, np.float32(2.0), 0.4, np.array([0.0625, 0.25, 0.375, 0.25, 0.0625], np.float32), *WORKING_LAB ), @@ -140,6 +143,7 @@ def _filled(fill): "sensor_log_unmix": "f7eb8efe1bbcfe884fd2e849e99d0546fac2f4829d3fb558d54a0128e681e53b", "lut_f32": "467bb9fd495d4797d912f6aeabd34f9a61b56491105464a0c852ed524488b1c1", "matrix_trc_u8": "acec1eadffb13dfef8473a02fb250b4ad8f57ac211a6265661b77d3cad49ab16", + "matrix_trc_u16": "8fd1af40718fc9bc59541e1128e46cb32e08ea05626f7cbd4cebe7905cecd02a", "usm_sharpen": "7267d2a3976c9381eadde01fe11ea60d46345b4747283e6ac68a2aa4a025d0a8", "sabattier": "f8956f67786c7c1d0e628af4b82fb12b6dde42e70da2504585f5bf20e7fbb461", "lith": "98dc05c62e191da2ad00c1f2327c48beff5af7e882be64cac5302cb0ea5099e5", diff --git a/tests/test_export_color_management.py b/tests/test_export_color_management.py index 22061749d..6174850ba 100644 --- a/tests/test_export_color_management.py +++ b/tests/test_export_color_management.py @@ -101,9 +101,13 @@ def test_cms_codec_unavailable_falls_back_to_lut_transform(self): # (hasattr/getattr on the module) re-enter that same __getattr__ and recurse. imagecodecs.cms_transform # noqa: B018 - with patch( - "negpy.services.rendering.image_processor.imagecodecs.cms_transform", - side_effect=ImportError("could not import name 'cms_transform' from 'imagecodecs'"), + # The bundled spaces are matrix/TRC and skip lcms; a LUT profile pair is what reaches it. + with ( + patch("negpy.services.rendering.image_processor.apply_matrix_trc_u16", return_value=None), + patch( + "negpy.services.rendering.image_processor.imagecodecs.cms_transform", + side_effect=ImportError("could not import name 'cms_transform' from 'imagecodecs'"), + ), ): out_fallback, icc = self.proc._apply_color_management_u16(img_u16, WORKING_COLOR_SPACE, ColorSpace.SRGB.value, None, None) @@ -132,9 +136,12 @@ def test_a_failed_transform_still_tags_the_pixels_it_ships(self): working_icc = icc_bytes_for_space(WORKING_COLOR_SPACE) imagecodecs.cms_transform # noqa: B018 (materialize before patching; see the test above) - with patch( - "negpy.services.rendering.image_processor.imagecodecs.cms_transform", - side_effect=RuntimeError("lcms2 refused the profile"), + with ( + patch("negpy.services.rendering.image_processor.apply_matrix_trc_u16", return_value=None), + patch( + "negpy.services.rendering.image_processor.imagecodecs.cms_transform", + side_effect=RuntimeError("lcms2 refused the profile"), + ), ): out, icc = self.proc._apply_color_management_u16(img_u16, WORKING_COLOR_SPACE, ColorSpace.SRGB.value, None, None) np.testing.assert_array_equal(out, img_u16) @@ -148,6 +155,21 @@ def test_a_failed_transform_still_tags_the_pixels_it_ships(self): np.testing.assert_array_equal(out, img_u16) self.assertEqual(icc, working_icc) + def test_16bit_matrix_trc_export_matches_exact_lcms(self): + """Every bundled export space is matrix/TRC: the 16-bit export takes the exact kernel, within + one level of lcms evaluated without its grid optimization, gamut edge included.""" + import imagecodecs + + img_u16 = np.random.default_rng(5).integers(0, 65536, size=(64, 64, 3), dtype=np.uint16) + src = icc_bytes_for_space(WORKING_COLOR_SPACE) + for target in (ColorSpace.SRGB.value, ColorSpace.P3_D65.value, ColorSpace.PROPHOTO.value, ColorSpace.REC2020.value): + with self.subTest(target=target): + out, icc = self.proc._apply_color_management_u16(img_u16, WORKING_COLOR_SPACE, target, None, None) + dst = icc_bytes_for_space(target) + exact = imagecodecs.cms_transform(img_u16, src, dst, colorspace="RGB", outcolorspace="RGB", intent=1, flags=0x2000 | 0x0100) + self.assertEqual(icc, dst) + self.assertLessEqual(int(np.abs(out.astype(np.int64) - exact.astype(np.int64)).max()), 1) + class TestDisplayTransform(unittest.TestCase): def test_srgb_working_is_identity(self): From 27f5f5d2fd130e20539a0d94632c670acb32c31c Mon Sep 17 00:00:00 2001 From: Marcin Zawalski Date: Sun, 25 Oct 2026 09:22:47 +0100 Subject: [PATCH 05/12] perf(export): greyscale export through a code table The greyscale export re-encoded every pixel with NumPy (working TRC decode, sRGB OETF, 16-bit, 8-bit), several full-size temporaries and np.power on each. The result depends only on the uint16 luma code, so the same formula now fills a 65536-entry table once per working space and negpy_core.lut_u16 reads it in parallel. The 8-bit path also drops its PIL round trip. Output is bit-identical (test against the per-pixel formula). 23.5 MP: 16-bit greyscale 207 -> 12 ms, 8-bit greyscale 246 -> 8.5 ms. --- negpy/services/rendering/image_processor.py | 39 ++++++++++++++------- negpy_core/src/convert.rs | 6 ++++ negpy_core/src/py/convert.rs | 30 ++++++++++++++++ tests/test_export_color_management.py | 22 ++++++++++++ 4 files changed, 84 insertions(+), 13 deletions(-) diff --git a/negpy/services/rendering/image_processor.py b/negpy/services/rendering/image_processor.py index 9a05f66cc..cdbb667c0 100644 --- a/negpy/services/rendering/image_processor.py +++ b/negpy/services/rendering/image_processor.py @@ -161,6 +161,23 @@ def run(i: int) -> None: return out +@lru_cache(maxsize=8) +def _grey_reencode_lut(working_color_space: Optional[str], bits: int) -> np.ndarray: + """The greyscale export re-encode as a table over every uint16 luma code: working TRC + decode, sRGB OETF, 16-bit, then 8-bit for ``bits`` 8. None keeps the working TRC.""" + codes = np.arange(65536, dtype=np.uint16) + out = codes + if working_color_space is not None: + lin = np.clip(np.asarray(working_oetf_decode(codes.astype(np.float32) / 65535.0)), 0.0, 1.0) + gray = np.where(lin <= 0.0031308, lin * 12.92, 1.055 * np.power(lin, 1.0 / 2.4) - 0.055) + out = np.clip(gray * 65535.0 + 0.5, 0.0, 65535.0).astype(np.uint16) + if bits == 8: + out = np.clip(np.round(out.astype(np.float32) / 257.0), 0.0, 255.0).astype(np.uint8) + out = np.ascontiguousarray(out) + out.setflags(write=False) + return out + + @lru_cache(maxsize=16) def _read_icc_bytes(path: str) -> bytes: with open(path, "rb") as f: @@ -1852,8 +1869,7 @@ def _export_pixels( return self._apply_color_management_u16_greyscale(img_int, working_color_space, color_space, icc_output, icc_input) return self._apply_color_management_u16(float_to_uint16(buffer), working_color_space, color_space, icc_output, icc_input) if is_greyscale: - pil_img, icc_bytes = self._greyscale_to_pil_u8(buffer, working_color_space, color_space, icc_output, icc_input) - return np.asarray(pil_img), icc_bytes + return self._greyscale_to_u8(buffer, working_color_space, color_space, icc_output, icc_input) return self.apply_color_management(float_to_uint8(buffer), working_color_space, color_space, icc_output, icc_input) def render_display_array( @@ -2132,30 +2148,27 @@ def _apply_color_management_u16_greyscale( """ if working_color_space == color_space: return img_u16, self._get_target_icc_bytes(color_space, output_icc_path) - lin = np.clip(np.asarray(working_oetf_decode(img_u16.astype(np.float32) / 65535.0)), 0.0, 1.0) - gray = np.where(lin <= 0.0031308, lin * 12.92, 1.055 * np.power(lin, 1.0 / 2.4) - 0.055) - out = np.clip(gray * 65535.0 + 0.5, 0.0, 65535.0).astype(np.uint16) + out = negpy_core.lut_u16(np.ascontiguousarray(img_u16), _grey_reencode_lut(working_color_space, 16)) return out, self._get_target_icc_bytes(color_space, output_icc_path) - def _greyscale_to_pil_u8( + def _greyscale_to_u8( self, buffer: np.ndarray, working_color_space: str, color_space: str, output_icc_path: Optional[str], input_icc_path: Optional[str] = None, - ) -> Tuple[Image.Image, Optional[bytes]]: - """Greyscale buffer -> color-managed 8-bit L image for JPEG/WebP. + ) -> Tuple[np.ndarray, Optional[bytes]]: + """Greyscale buffer -> color-managed 8-bit luma for JPEG/WebP. Runs the 16-bit grey re-encode first so the 8-bit result carries the same TRC (and profile bytes) as the greyscale TIFF/PNG of the same edit. """ img16 = float_to_uint_luma(np.ascontiguousarray(buffer), bit_depth=16) - img16, icc_bytes = self._apply_color_management_u16_greyscale( - img16, working_color_space, color_space, output_icc_path, input_icc_path - ) - img8 = np.clip(np.round(img16.astype(np.float32) / 257.0), 0.0, 255.0).astype(np.uint8) - return Image.fromarray(img8), icc_bytes + icc_bytes = self._get_target_icc_bytes(color_space, output_icc_path) + if working_color_space == color_space: + return negpy_core.lut_u16(img16, _grey_reencode_lut(None, 8)), icc_bytes + return negpy_core.lut_u16(img16, _grey_reencode_lut(working_color_space, 8)), icc_bytes def _apply_color_management_u16( self, diff --git a/negpy_core/src/convert.rs b/negpy_core/src/convert.rs index 120f0ad6b..ad8a9c32e 100644 --- a/negpy_core/src/convert.rs +++ b/negpy_core/src/convert.rs @@ -47,6 +47,12 @@ pub fn scaled_to_float32 + Sync>(src: &[T], scale: f32, dst: crate::par_zip(dst, src, move |d, &v| *d = v.into() * scale); } +/// `lut[src]` of a uint16 array through a 65536-entry table. +pub fn lut_u16(src: &[u16], lut: &[T], dst: &mut [T]) { + let lut = &lut[..65536]; + crate::par_zip(dst, src, move |d, &v| *d = lut[v as usize]); +} + /// Numba's fastmath luminance in float64: two fused multiply-adds. fn luma(p: &[f32]) -> f64 { LUMA_B.mul_add( diff --git a/negpy_core/src/py/convert.rs b/negpy_core/src/py/convert.rs index 58b33a68c..a061f7311 100644 --- a/negpy_core/src/py/convert.rs +++ b/negpy_core/src/py/convert.rs @@ -360,6 +360,35 @@ fn int_to_unit_float32<'py>( .into_any()) } +/// `lut[img]` of a contiguous uint16 array through a 65536-entry uint8 or uint16 table. +#[pyfunction] +fn lut_u16<'py>( + py: Python<'py>, + img: PyReadonlyArrayDyn<'py, u16>, + lut: &Bound<'py, PyAny>, +) -> PyResult> { + let src = contiguous(&img)?; + if let Ok(t) = lut.extract::>() { + let table = contiguous(&t)?; + if table.len() != 65536 { + return Err(PyValueError::new_err("needs a 65536-entry table")); + } + return Ok(filled(py, img.shape(), "uint16", |d| { + convert::lut_u16(src, table, d) + })? + .into_any()); + } + let t = lut.extract::>()?; + let table = contiguous(&t)?; + if table.len() != 65536 { + return Err(PyValueError::new_err("needs a 65536-entry table")); + } + Ok(filled(py, img.shape(), "uint8", |d| { + convert::lut_u16(src, table, d) + })? + .into_any()) +} + /// `img`, an HxWxP uint8 view of an image with P-byte pixels, turned by EXIF `orientation`. #[pyfunction] fn orient<'py>( @@ -382,6 +411,7 @@ fn orient<'py>( pub fn register(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_function(wrap_pyfunction!(orient, m)?)?; + m.add_function(wrap_pyfunction!(lut_u16, m)?)?; m.add_function(wrap_pyfunction!(int_to_unit_float32, m)?)?; m.add_function(wrap_pyfunction!(nan_to_num, m)?)?; m.add_function(wrap_pyfunction!(gray_to_rgb, m)?)?; diff --git a/tests/test_export_color_management.py b/tests/test_export_color_management.py index 6174850ba..24f8af85a 100644 --- a/tests/test_export_color_management.py +++ b/tests/test_export_color_management.py @@ -55,6 +55,28 @@ def test_greyscale_export_matches_gamma22_tag(self): decoded = (out.astype(np.float32) / 65535.0) ** 2.2 # as a Gamma2.2 viewer would decode np.testing.assert_allclose(decoded.reshape(-1), lin, atol=0.01) + def test_greyscale_reencode_table_matches_the_per_pixel_formula(self): + """The 16- and 8-bit greyscale exports read one table per code; it must give the bits the + per-pixel re-encode gives.""" + from negpy.kernel.image.logic import working_oetf_decode + + u16 = np.random.default_rng(9).integers(0, 65536, size=(97, 131), dtype=np.uint16) + lin = np.clip(np.asarray(working_oetf_decode(u16.astype(np.float32) / 65535.0)), 0.0, 1.0) + gray = np.where(lin <= 0.0031308, lin * 12.92, 1.055 * np.power(lin, 1.0 / 2.4) - 0.055) + expected16 = np.clip(gray * 65535.0 + 0.5, 0.0, 65535.0).astype(np.uint16) + out16, _ = self.proc._apply_color_management_u16_greyscale(u16, WORKING_COLOR_SPACE, ColorSpace.GREYSCALE.value, None) + np.testing.assert_array_equal(out16, expected16) + + buffer = np.random.default_rng(10).random((61, 83), dtype=np.float32) + from negpy.kernel.image.logic import float_to_uint_luma + + img16 = float_to_uint_luma(buffer, bit_depth=16) + for target in (ColorSpace.GREYSCALE.value, WORKING_COLOR_SPACE): + mid, _ = self.proc._apply_color_management_u16_greyscale(img16, WORKING_COLOR_SPACE, target, None) + expected8 = np.clip(np.round(mid.astype(np.float32) / 257.0), 0.0, 255.0).astype(np.uint8) + out8, _ = self.proc._greyscale_to_u8(buffer, WORKING_COLOR_SPACE, target, None) + np.testing.assert_array_equal(out8, expected8) + def test_export_is_appearance_preserving_across_spaces(self): """Exporting to sRGB / Adobe / ProPhoto must look the same in a CM viewer. From 73dc4807db0f428ffd6e5e142c60505778490853 Mon Sep 17 00:00:00 2001 From: Marcin Zawalski Date: Sun, 25 Oct 2026 09:27:01 +0100 Subject: [PATCH 06/12] perf(core): Burn vignette in Rust apply_vignette was the last NumPy pixel stage in the CPU finish: a full meshgrid, sqrt, cos and exp2 with about ten full-size temporaries. The CPU engine renders every thumbnail refresh, so it ran there as well as in the CPU preview and export. negpy_core.apply_vignette keeps NumPy's float32 order with cos and exp2 from libm in float64, rounded once: within 4 ULP of NumPy's SVML results, the same bits on every platform (same-bits pin). THUMBNAIL_RENDER_VERSION 17 -> 18 for those ULPs. 24 threads: 1600 px 31 -> 6 ms, 23.5 MP 386 -> 77 ms. --- negpy/features/finish/logic.py | 24 +--------- .../services/assets/thumbnail_fingerprint.py | 2 +- negpy_core/src/finish.rs | 45 +++++++++++++++++++ negpy_core/src/py/finish.rs | 29 ++++++++++++ tests/test_core_same_bits.py | 2 + tests/thumbnail_fingerprint_defaults.json | 2 +- 6 files changed, 79 insertions(+), 25 deletions(-) diff --git a/negpy/features/finish/logic.py b/negpy/features/finish/logic.py index 9ba36ce16..71d8a7f14 100644 --- a/negpy/features/finish/logic.py +++ b/negpy/features/finish/logic.py @@ -114,26 +114,4 @@ def apply_vignette(img: ImageBuffer, stops: float, size: float, roundness: float """ if stops == 0.0: return img - - h, w = img.shape[:2] - cy, cx = (h - 1) * 0.5, (w - 1) * 0.5 - - y_coords = np.arange(h, dtype=np.float32) - x_coords = np.arange(w, dtype=np.float32) - yy, xx = np.meshgrid(y_coords, x_coords, indexing="ij") - max_dist = float(np.sqrt(cx * cx + cy * cy)) - d_radial = np.sqrt((xx - cx) ** 2 + (yy - cy) ** 2) / max(max_dist, 1.0) - d_rect = np.maximum(np.abs(xx - cx) / max(cx, 1.0), np.abs(yy - cy) / max(cy, 1.0)) - dist = d_radial * (1.0 - roundness) + d_rect * roundness - - # Remap: size=0 → vignette barely at corners, size=1 → covers entire image - midpoint = 1.0 - size - t = (dist - midpoint) / max(1.0 - midpoint, 1e-6) - t = np.clip(t, 0.0, 1.0) - - # Smooth cosine falloff - factor = 0.5 * (1.0 - np.cos(t * np.pi)) - - result = img * np.exp2(-stops * factor[:, :, np.newaxis]) - - return ensure_image(np.clip(result, 0.0, 1.0)) + return ensure_image(negpy_core.apply_vignette(np.ascontiguousarray(img, dtype=np.float32), float(stops), float(size), float(roundness))) diff --git a/negpy/services/assets/thumbnail_fingerprint.py b/negpy/services/assets/thumbnail_fingerprint.py index e2566ab18..3e6048334 100644 --- a/negpy/services/assets/thumbnail_fingerprint.py +++ b/negpy/services/assets/thumbnail_fingerprint.py @@ -12,7 +12,7 @@ from negpy.domain.models import WorkspaceConfig # Bump in any change that alters rendered pixels for unchanged settings, a changed default included. -THUMBNAIL_RENDER_VERSION = 17 +THUMBNAIL_RENDER_VERSION = 18 # Marks a thumbnail made from the source preview, which runs none of the frame's settings. QUICK = "quick" diff --git a/negpy_core/src/finish.rs b/negpy_core/src/finish.rs index 33e27c5c7..27ba158c5 100644 --- a/negpy_core/src/finish.rs +++ b/negpy_core/src/finish.rs @@ -344,3 +344,48 @@ pub fn apply_carrier(src: &[f32], (h, w): (usize, usize), c: &Carrier, dst: &mut } }); } + +/// `apply_vignette`: an exposure change of `-stops * factor` stops, the factor a cosine ramp +/// over a radial/rectangular distance from the centre, in NumPy's float32 order. cos and exp2 +/// come from libm in float64, rounded once. +#[allow(clippy::too_many_arguments)] +pub fn apply_vignette( + src: &[f32], + h: usize, + w: usize, + cn: usize, + stops: f64, + size: f64, + roundness: f64, + dst: &mut [f32], +) { + use rayon::prelude::*; + let (cy, cx) = ((h as f64 - 1.0) * 0.5, (w as f64 - 1.0) * 0.5); + let max_dist = (cx * cx + cy * cy).sqrt().max(1.0) as f32; + let (cxf, cyf) = (cx as f32, cy as f32); + let (rx, ry) = (cx.max(1.0) as f32, cy.max(1.0) as f32); + let (keep, mix) = ((1.0 - roundness) as f32, roundness as f32); + let midpoint = 1.0 - size; + let (mid, span) = (midpoint as f32, (1.0 - midpoint).max(1e-6) as f32); + let neg_stops = -stops as f32; + let pi = std::f64::consts::PI as f32; + dst.par_chunks_mut(w * cn) + .zip(src.par_chunks(w * cn)) + .enumerate() + .for_each(move |(y, (d, s))| { + let dy = y as f32 - cyf; + for x in 0..w { + let dx = x as f32 - cxf; + let d_radial = (dx * dx + dy * dy).sqrt() / max_dist; + let d_rect = (dx.abs() / rx).max(dy.abs() / ry); + let dist = d_radial * keep + d_rect * mix; + let t = ((dist - mid) / span).clamp(0.0, 1.0); + let factor = 0.5f32 * (1.0f32 - libm::cos((t * pi) as f64) as f32); + let gain = libm::exp2((neg_stops * factor) as f64) as f32; + for c in 0..cn { + let i = x * cn + c; + d[i] = (s[i] * gain).clamp(0.0, 1.0); + } + } + }); +} diff --git a/negpy_core/src/py/finish.rs b/negpy_core/src/py/finish.rs index c3ee0ab9d..e91630e20 100644 --- a/negpy_core/src/py/finish.rs +++ b/negpy_core/src/py/finish.rs @@ -54,6 +54,35 @@ fn apply_carrier<'py>( Ok(out) } +/// `apply_vignette(img, stops, size, roundness)` of a contiguous HxWxC float32 image. +#[pyfunction] +fn apply_vignette<'py>( + py: Python<'py>, + img: PyReadonlyArrayDyn<'py, f32>, + stops: f64, + size: f64, + roundness: f64, +) -> PyResult>> { + let [h, w, cn] = img.shape()[..] else { + return Err(PyValueError::new_err("needs an HxWxC image")); + }; + let src = img + .as_slice() + .map_err(|_| PyValueError::new_err("needs a contiguous image"))?; + let out = PyArrayDyn::::zeros(py, img.shape(), false); + { + let mut rw = out.readwrite(); + let dst = rw + .as_slice_mut() + .map_err(|e| PyValueError::new_err(e.to_string()))?; + super::detach(py, dst.len(), || { + finish::apply_vignette(src, h, w, cn, stops, size, roundness, dst) + }); + } + Ok(out) +} + pub fn register(m: &Bound<'_, PyModule>) -> PyResult<()> { + m.add_function(wrap_pyfunction!(apply_vignette, m)?)?; m.add_function(wrap_pyfunction!(apply_carrier, m)?) } diff --git a/tests/test_core_same_bits.py b/tests/test_core_same_bits.py index 3a94471d0..c1f9767af 100644 --- a/tests/test_core_same_bits.py +++ b/tests/test_core_same_bits.py @@ -73,6 +73,7 @@ def _filled(fill): ), "lut_f32": lambda: negpy_core.apply_lut_f32(_IMG3, _LUT), "matrix_trc_u8": lambda: negpy_core.matrix_trc_u8((_IMG3 * 255).astype(np.uint8), np.ascontiguousarray(_RAMP[:, ::256]), _M, _RAMP), + "vignette": lambda: negpy_core.apply_vignette(_IMG3, 0.7, 0.6, 0.3), "matrix_trc_u16": lambda: negpy_core.matrix_trc_u16( (_IMG3 * 65535).astype(np.uint16), np.ascontiguousarray(np.tile(_LIN, (3, 1))), _M, np.ascontiguousarray(_RAMP[:, ::16]) ), @@ -144,6 +145,7 @@ def _filled(fill): "lut_f32": "467bb9fd495d4797d912f6aeabd34f9a61b56491105464a0c852ed524488b1c1", "matrix_trc_u8": "acec1eadffb13dfef8473a02fb250b4ad8f57ac211a6265661b77d3cad49ab16", "matrix_trc_u16": "8fd1af40718fc9bc59541e1128e46cb32e08ea05626f7cbd4cebe7905cecd02a", + "vignette": "dea377797d6b0852a8dec1665c95aa87c9ab2b6a8e03f808e27140bc4b6566c6", "usm_sharpen": "7267d2a3976c9381eadde01fe11ea60d46345b4747283e6ac68a2aa4a025d0a8", "sabattier": "f8956f67786c7c1d0e628af4b82fb12b6dde42e70da2504585f5bf20e7fbb461", "lith": "98dc05c62e191da2ad00c1f2327c48beff5af7e882be64cac5302cb0ea5099e5", diff --git a/tests/thumbnail_fingerprint_defaults.json b/tests/thumbnail_fingerprint_defaults.json index 9e86d3e03..768b3ea84 100644 --- a/tests/thumbnail_fingerprint_defaults.json +++ b/tests/thumbnail_fingerprint_defaults.json @@ -206,5 +206,5 @@ "vanadium_strength": "0.0" } }, - "render_version": 17 + "render_version": 18 } From 870e7b9d32930afe975da493356edfc502de790e Mon Sep 17 00:00:00 2001 From: Marcin Zawalski Date: Sun, 25 Oct 2026 09:30:40 +0100 Subject: [PATCH 07/12] perf(core): flat master curve in Rust apply_flat_curve, the one pixel step of the Flat intent, ran in NumPy: an offset add, a BLAS product of the full image with the stock matrix and the affine code map, each a full-size temporary. The flat export is CPU-only. negpy_core.flat_curve does it in one parallel pass, the matrix in matmul_3x3's order; the output is bit-identical to the NumPy version here. 23.5 MP: 409 -> 24 ms per flat export. --- .github/workflows/ci.yml | 2 +- negpy/features/exposure/logic.py | 8 +++--- negpy_core/src/color.rs | 43 ++++++++++++++++++++++++++++++++ negpy_core/src/py/color.rs | 28 +++++++++++++++++++++ tests/test_core_flat_curve.py | 29 +++++++++++++++++++++ tests/test_core_same_bits.py | 2 ++ 6 files changed, 106 insertions(+), 6 deletions(-) create mode 100644 tests/test_core_flat_curve.py diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 799079d61..600acce41 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -129,4 +129,4 @@ jobs: run: uv sync --group dev --group plustek - name: negpy_core pins and oracles - run: uv run pytest tests/test_percentiles.py tests/test_core_same_bits.py tests/test_core_orient.py tests/test_core_ir_downsample.py tests/test_core_resize.py tests/test_core_filter.py tests/test_core_morph.py tests/test_core_gray.py tests/test_core_contours.py tests/test_core_shapes.py tests/test_core_components.py tests/test_core_color.py tests/test_core_warp.py tests/test_core_lab.py tests/test_core_altprocess.py + run: uv run pytest tests/test_percentiles.py tests/test_core_same_bits.py tests/test_core_orient.py tests/test_core_ir_downsample.py tests/test_core_flat_curve.py tests/test_core_resize.py tests/test_core_filter.py tests/test_core_morph.py tests/test_core_gray.py tests/test_core_contours.py tests/test_core_shapes.py tests/test_core_components.py tests/test_core_color.py tests/test_core_warp.py tests/test_core_lab.py tests/test_core_altprocess.py diff --git a/negpy/features/exposure/logic.py b/negpy/features/exposure/logic.py index 7cb0e151a..b93d602cf 100644 --- a/negpy/features/exposure/logic.py +++ b/negpy/features/exposure/logic.py @@ -558,11 +558,9 @@ def apply_flat_curve( editing. WB rides as an additive per-channel shift in log space; the print stock's density matrix (rows sum 1) mixes after it, as in the print kernel. """ - val = np.asarray(image, dtype=np.float32) + np.asarray(cmy_offsets, dtype=np.float32) - if stock_matrix is not None: - val = val @ np.asarray(stock_matrix, dtype=np.float32).T - code = lift + gain * (1.0 - val) - return ensure_image(np.clip(code, 0.0, 1.0).astype(np.float32)) + m = None if stock_matrix is None else np.ascontiguousarray(stock_matrix, dtype=np.float32) + offsets = tuple(float(np.float32(o)) for o in cmy_offsets) + return ensure_image(negpy_core.flat_curve(np.ascontiguousarray(image, dtype=np.float32), float(gain), float(lift), offsets, m)) def default_grade_range() -> float: diff --git a/negpy_core/src/color.rs b/negpy_core/src/color.rs index 596fa9544..1e408a652 100644 --- a/negpy_core/src/color.rs +++ b/negpy_core/src/color.rs @@ -337,6 +337,49 @@ fn matmul_block(src: &[f32], dst: &mut [f32], m: [[f32; 3]; 3]) { } } +/// `apply_flat_curve`: `lift + gain * (1 - v)` clipped to 0..1, with v the pixel plus the +/// per-channel `offsets`, mixed by the 3x3 `m` (rows sum to 1) in `matmul_3x3`'s order when given. +pub fn flat_curve( + src: &[f32], + dst: &mut [f32], + offsets: [f32; 3], + m: Option<[[f32; 3]; 3]>, + gain: f32, + lift: f32, +) { + dst.par_chunks_mut(BLOCK * 3) + .zip(src.par_chunks(BLOCK * 3)) + .for_each(|(d, p)| flat_curve_block(p, d, offsets, m, gain, lift)); +} + +#[multiversion(targets( + "x86_64+avx512f+avx512bw+avx512vl+avx512dq", + "x86_64+avx2+fma", + "aarch64+neon" +))] +fn flat_curve_block( + src: &[f32], + dst: &mut [f32], + offsets: [f32; 3], + m: Option<[[f32; 3]; 3]>, + gain: f32, + lift: f32, +) { + for (d, p) in dst + .as_chunks_mut::<3>() + .0 + .iter_mut() + .zip(src.as_chunks::<3>().0) + { + let [r, g, b] = std::array::from_fn(|c| p[c] + offsets[c]); + let v = match m { + Some(m) => m.map(|row| b.mul_add(row[2], g.mul_add(row[1], row[0] * r))), + None => [r, g, b], + }; + *d = v.map(|x| (lift + gain * (1.0 - x)).clamp(0.0, 1.0)); + } +} + /// `_in_gamut_lab`: (L, a, b) decodes to linear RGB within 0..1, 1e-4 either side; the bounds /// are float32 sums, compared in float64. #[inline(always)] diff --git a/negpy_core/src/py/color.rs b/negpy_core/src/py/color.rs index 9455309e8..ea79edb32 100644 --- a/negpy_core/src/py/color.rs +++ b/negpy_core/src/py/color.rs @@ -166,7 +166,35 @@ fn gamut_aware_chroma_scale<'py>( }) } +/// `apply_flat_curve(img, gain, lift, offsets, stock_matrix)` of an (..., 3) float32 array. +#[pyfunction] +#[pyo3(signature = (px, gain, lift, offsets, m=None))] +fn flat_curve<'py>( + py: Python<'py>, + px: PyReadonlyArrayDyn<'py, f32>, + gain: f32, + lift: f32, + offsets: [f32; 3], + m: Option>, +) -> PyResult>> { + triplets(&px)?; + let m = match m { + Some(m) => { + let m = m.as_array(); + if m.shape() != [3, 3] { + return Err(PyValueError::new_err("needs a 3x3 matrix")); + } + Some(std::array::from_fn(|r| std::array::from_fn(|c| m[[r, c]]))) + } + None => None, + }; + filled(py, &px, |src, d| { + color::flat_curve(src, d, offsets, m, gain, lift) + }) +} + pub fn register(m: &Bound<'_, PyModule>) -> PyResult<()> { + m.add_function(wrap_pyfunction!(flat_curve, m)?)?; m.add_function(wrap_pyfunction!(oetf_encode, m)?)?; m.add_function(wrap_pyfunction!(oetf_decode, m)?)?; m.add_function(wrap_pyfunction!(rgb_to_lab, m)?)?; diff --git a/tests/test_core_flat_curve.py b/tests/test_core_flat_curve.py new file mode 100644 index 000000000..def26df97 --- /dev/null +++ b/tests/test_core_flat_curve.py @@ -0,0 +1,29 @@ +import numpy as np + +from negpy.features.exposure.logic import apply_flat_curve + +_IMG = np.random.default_rng(4).random((61, 89, 3), dtype=np.float32) * 1.4 - 0.2 +_M = np.array([[0.8, 0.15, 0.05], [0.1, 0.85, 0.05], [0.02, 0.08, 0.9]], np.float32) + + +def _numpy(img, gain, lift, offsets, m): + val = img + np.asarray(offsets, dtype=np.float32) + if m is not None: + val = val @ m.T + return np.clip(lift + gain * (1.0 - val), 0.0, 1.0).astype(np.float32) + + +def test_flat_curve_without_a_stock_is_numpys_arithmetic(): + out = apply_flat_curve(_IMG, 0.65, 0.10, (0.013, -0.02, 0.0)) + assert np.array_equal(out, _numpy(_IMG, 0.65, 0.10, (0.013, -0.02, 0.0), None)) + + +def test_flat_curve_with_a_stock_matches_the_matrix_product(): + out = apply_flat_curve(_IMG, 0.65, 0.10, (0.0, 0.01, -0.01), _M) + np.testing.assert_allclose(out, _numpy(_IMG, 0.65, 0.10, (0.0, 0.01, -0.01), _M), atol=2e-7) + + +def test_flat_curve_keeps_nan_and_clips(): + img = np.array([[[np.nan, -5.0, 5.0]]], np.float32) + out = apply_flat_curve(img, 0.65, 0.10) + assert np.isnan(out[0, 0, 0]) and out[0, 0, 1] == 1.0 and out[0, 0, 2] == 0.0 diff --git a/tests/test_core_same_bits.py b/tests/test_core_same_bits.py index c1f9767af..4176a62e1 100644 --- a/tests/test_core_same_bits.py +++ b/tests/test_core_same_bits.py @@ -74,6 +74,7 @@ def _filled(fill): "lut_f32": lambda: negpy_core.apply_lut_f32(_IMG3, _LUT), "matrix_trc_u8": lambda: negpy_core.matrix_trc_u8((_IMG3 * 255).astype(np.uint8), np.ascontiguousarray(_RAMP[:, ::256]), _M, _RAMP), "vignette": lambda: negpy_core.apply_vignette(_IMG3, 0.7, 0.6, 0.3), + "flat_curve": lambda: negpy_core.flat_curve(_IMG3, 0.65, 0.1, (0.013, -0.02, 0.0), _M), "matrix_trc_u16": lambda: negpy_core.matrix_trc_u16( (_IMG3 * 65535).astype(np.uint16), np.ascontiguousarray(np.tile(_LIN, (3, 1))), _M, np.ascontiguousarray(_RAMP[:, ::16]) ), @@ -146,6 +147,7 @@ def _filled(fill): "matrix_trc_u8": "acec1eadffb13dfef8473a02fb250b4ad8f57ac211a6265661b77d3cad49ab16", "matrix_trc_u16": "8fd1af40718fc9bc59541e1128e46cb32e08ea05626f7cbd4cebe7905cecd02a", "vignette": "dea377797d6b0852a8dec1665c95aa87c9ab2b6a8e03f808e27140bc4b6566c6", + "flat_curve": "e47f6084a58e19119c6e3364ebf017dfc7f039ec6c73e0eea8e942519021bf0b", "usm_sharpen": "7267d2a3976c9381eadde01fe11ea60d46345b4747283e6ac68a2aa4a025d0a8", "sabattier": "f8956f67786c7c1d0e628af4b82fb12b6dde42e70da2504585f5bf20e7fbb461", "lith": "98dc05c62e191da2ad00c1f2327c48beff5af7e882be64cac5302cb0ea5099e5", From fa244408fad2b8f6480c193adba7d3df2a400192 Mon Sep 17 00:00:00 2001 From: Marcin Zawalski Date: Sun, 25 Oct 2026 09:46:19 +0100 Subject: [PATCH 08/12] perf(core): slide transfer curve in Rust apply_transfer_curve, the CPU render of every Slide and Positive frame, was the one render stage the port missed: NumPy with dozens of full-size temporaries, softplus, tanh and power per channel. It runs for every slide thumbnail refresh (always CPU), the CPU preview and CPU exports. negpy_core.apply_transfer_curve keeps NumPy's float32 order and runs each step as a flat loop over runs of 256 pixels, so the mathf float64 transcendentals vectorize; they round once to float32, within 3e-7 of NumPy's SVML results and the same bits on every platform (same-bits pin). The NumPy helpers it replaced (_softplus, _sigmoid, _black_taper, separation_damping_gain_np) are gone. Covered by the existing render version bump on this branch. 24 threads, 1600 px: defaults 57 -> 3.5 ms, toe/shoulder/contrast 251 -> 8.5 ms, with Dye Separation 725 -> 20 ms. 23.5 MP edited: 4.0 s -> 0.11 s. --- negpy/features/exposure/logic.py | 11 - negpy/features/transparency/logic.py | 164 ++++-------- .../transparency/shaders/transfer.wgsl | 2 +- negpy_core/src/lib.rs | 2 + negpy_core/src/py/mod.rs | 1 + negpy_core/src/py/transfer.rs | 98 +++++++ negpy_core/src/transfer.rs | 240 ++++++++++++++++++ tests/test_core_same_bits.py | 15 ++ 8 files changed, 403 insertions(+), 130 deletions(-) create mode 100644 negpy_core/src/py/transfer.rs create mode 100644 negpy_core/src/transfer.rs diff --git a/negpy/features/exposure/logic.py b/negpy/features/exposure/logic.py index b93d602cf..28f89c84c 100644 --- a/negpy/features/exposure/logic.py +++ b/negpy/features/exposure/logic.py @@ -59,17 +59,6 @@ def separation_damping_gain(k: float, damping: float, chroma: float, ref_spread: return float(kf) -def separation_damping_gain_np(k: float, damping: float, chroma: Any, ref_spread: float) -> Any: - """Vectorized numpy twin of separation_damping_gain, for the transfer curve's - per-pixel gain over a whole image rather than the print curve's per-pixel scalar. - Same math; see that function for the derivation.""" - if k <= 0.0: - return np.zeros_like(chroma, dtype=np.float32) - h = (np.float32(ref_spread) - chroma) / (np.float32(ref_spread) + chroma) - kf = np.power(np.float32(k), (1.0 - np.float32(damping)) + np.float32(damping) * h, dtype=np.float32) - return np.minimum(kf, np.float32(3.0)) - - def tone_key_weight_np(lum: Any, e0: float, e1: float) -> Any: """A tone-limited mask's weight at luma `lum`, over a whole raster: a smoothstep from 0 at e0 to 1 at e1 (see placement.key_edges). The print curve in negpy_core and exposure.wgsl mirror it.""" diff --git a/negpy/features/transparency/logic.py b/negpy/features/transparency/logic.py index 9585fdea4..55edff22b 100644 --- a/negpy/features/transparency/logic.py +++ b/negpy/features/transparency/logic.py @@ -35,6 +35,7 @@ from typing import Any, Dict, Mapping, Optional, Tuple +import negpy_core import numpy as np from negpy.domain.types import ImageBuffer @@ -45,7 +46,6 @@ per_channel_dye_separation, per_channel_toe_shoulder, per_channel_widths, - separation_damping_gain_np, slope_to_grade, ) from negpy.features.exposure.models import EXPOSURE_CONSTANTS, ExposureConfig @@ -155,31 +155,9 @@ def display_rendering(scene_linear: np.ndarray) -> np.ndarray: return np.clip(num / np.maximum(den, np.float32(1e-8)), 0.0, 1.0).astype(np.float32) -def _sigmoid(x: np.ndarray) -> np.ndarray: - """Logistic sigmoid, overflow-safe for large |x| (mirrors logic.py::_fast_sigmoid).""" - return (0.5 * (1.0 + np.tanh(0.5 * np.asarray(x, dtype=np.float32)))).astype(np.float32) - - ZONE_BLACK_TAPER = 1.0 -def _black_taper(d: np.ndarray, density_range: float) -> np.ndarray: - """Fades a shadow lift back to nothing at the bottom of the window. - - On a print, Zone Density is bounded by paper black -- a shadow burn cannot exceed - d_max. This curve has no paper, so without a bound a lift walks the black point up - with it and the frame simply stops having blacks. Smoothstep so the taper adds no kink. - """ - t = np.clip((np.float32(density_range) - d) / np.float32(ZONE_BLACK_TAPER), 0.0, 1.0) - return (t * t * (np.float32(3.0) - np.float32(2.0) * t)).astype(np.float32) - - -def _softplus(x: np.ndarray, width: float) -> np.ndarray: - """width * log(1 + exp(x / width)), overflow-safe: ~x for x >> width, ~0 for x << -width.""" - t = x / width - return (width * (np.logaddexp(0.0, -np.abs(t)) + np.maximum(t, 0.0))).astype(np.float32) - - def transfer_widths(config: ExposureConfig) -> Tuple[Tuple[float, float, float], Tuple[float, float, float]]: """ Per-channel knee softness, in density, from the Toe/Shoulder Width sliders. Scaled @@ -424,102 +402,52 @@ def apply_transfer_curve( """ c = TRANSFER_CONSTANTS base_width = float(c["transfer_knee_width"]) - tw3 = toe_widths or (base_width, base_width, base_width) - sw3 = shoulder_widths or (base_width, base_width, base_width) - pivot = float(c["transfer_contrast_pivot"]) - toe_knee = float(c["transfer_toe_knee"]) - sh_knee = float(c["transfer_shoulder_knee"]) - - def shape(d: np.ndarray, ch: int) -> np.ndarray: - """One channel's capture density through every control above Dye Separation.""" - - # Cast Removal: the channel's neutral refs onto green's. First, so every control - # below shapes the corrected signal. - if cast_gain[ch] != 1.0 or cast_offset[ch] != 0.0: - d = d * np.float32(cast_gain[ch]) + np.float32(cast_offset[ch]) - - if exposure_offset != 0.0 or cmy_offsets[ch] != 0.0: - # cmy_offsets arrive in normalized space, because filtration_offsets divides by the - # channel's stretch range, so they scale back up by the same range. - d = d - np.float32(exposure_offset) + np.float32(cmy_offsets[ch] * density_range) - - if contrast != 1.0: - d = np.float32(pivot) + (d - np.float32(pivot)) * np.float32(contrast) - - # Shadows/Highlights WB: regional CMY, using the print path's own kernel - # (logic.py, "Regional CMY") and geometry (wb_split_geometry). - if shadow_cmy[ch] != 0.0 or highlight_cmy[ch] != 0.0: - wb_c, wb_k = wb_split_geometry() - w_sh = _sigmoid(np.float32(wb_k) * (d - np.float32(wb_c))) - w_hi = np.float32(1.0) - w_sh - d = d + np.float32(shadow_cmy[ch]) * w_sh + np.float32(highlight_cmy[ch]) * w_hi - - # Zone Density: mid-sparing brightness offsets, using the print path's own kernel and - # weights (logic.py, "Zone Density (ΔD)"). Positive adds density, so it darkens: the - # Density convention, the opposite sign to a Lightroom Shadows slider. Runs after - # contrast and before the knees, as it does on the print. - if shadow_density != 0.0 or highlight_density != 0.0: - sh_c, hi_c, k = zone_geometry() - w_sh = _sigmoid(np.float32(k) * (d - np.float32(sh_c))) - w_sh = w_sh * _black_taper(d, density_range) - w_hi = np.float32(1.0) - _sigmoid(np.float32(k) * (d - np.float32(hi_c))) - d = d + np.float32(shadow_density) * w_sh + np.float32(highlight_density) * w_hi - - # Shadows sit at high density and highlights at low, so the toe compresses above its knee - # and the shoulder below its own. - if toe[ch] != 0.0: - d = d - np.float32(toe[ch]) * _softplus(d - np.float32(toe_knee), tw3[ch]) - if shoulder[ch] != 0.0: - d = d + np.float32(shoulder[ch]) * _softplus(np.float32(sh_knee) - d, sw3[ch]) - - return d - - n = np.asarray(img_norm, dtype=np.float32) - capture = n * np.float32(density_range) - dens = np.empty_like(n) - for ch in range(3): - dens[:, :, ch] = shape(capture[:, :, ch], ch) - + wb_c, wb_k = wb_split_geometry() + sh_c, hi_c, zone_k = zone_geometry() sep_k3 = per_channel_dye_separation(separation, separation_trims) - if sep_k3 != (1.0, 1.0, 1.0): - # Both softplus steps rise with their input, so d_ref never falls as a channel - # darkens; a reference that can fall folds a smooth gradient into bands. The cap is - # decided on the capture, where the window clamp always sits at the dense edge, so - # no edit can move a clamped channel back under it. - d_lo = capture.min(axis=2, keepdims=True) - cap = np.float32(SEPARATION_CAP_LEVEL) + _softplus( - d_lo + np.float32(SEPARATION_CAP_SPREAD - SEPARATION_CAP_LEVEL), SEPARATION_CAP_SOFTNESS - ) - capped = cap - _softplus(cap - capture, SEPARATION_CAP_SOFTNESS) - capped = np.stack([shape(capped[:, :, ch], ch) for ch in range(3)], axis=2) - d_ref = capped.mean(axis=2, keepdims=True) - e = dens - d_ref - if damping > 0.0: - # Separation Damping makes each channel's k chroma-dependent per pixel, from - # the same chroma but each channel's own k (see separation_damping_gain_np). - # The chroma reads the capped channels too. - ce = capped - d_ref - chroma = np.sqrt(((ce[:, :, 0] - ce[:, :, 1]) ** 2 + (ce[:, :, 1] - ce[:, :, 2]) ** 2 + (ce[:, :, 0] - ce[:, :, 2]) ** 2) / 3.0) - ref_spread = float(EXPOSURE_CONSTANTS["separation_damping_ref_spread"]) - k_eff = np.stack( - [separation_damping_gain_np(sep_k3[ch], damping, chroma, ref_spread) for ch in range(3)], - axis=2, - ) - dens = d_ref + k_eff * e - else: - k3 = np.asarray(sep_k3, dtype=np.float32) - dens = d_ref + k3[np.newaxis, np.newaxis, :] * e - - out = np.power(np.float32(10.0), -dens, dtype=np.float32) - - # Baseline and display rendering last, so the controls above shape the scene and this - # only decides how the scene is shown. A finished positive sits nowhere below a sensor - # white level, so neither belongs: it passes through at scene scale, clipped to the - # same range display_rendering outputs. - if positive_source: - return ensure_image(np.clip(out, 0.0, 1.0)) - gain = np.float32(2.0 ** float(c["transfer_baseline_ev"])) - return ensure_image(display_rendering(out * gain)) + params = { + "density_range": density_range, + "cast_gain": tuple(cast_gain), + "cast_offset": tuple(cast_offset), + "cast_on": tuple(cast_gain[ch] != 1.0 or cast_offset[ch] != 0.0 for ch in range(3)), + "exposure_offset": exposure_offset, + # cmy_offsets arrive in normalized space, because filtration_offsets divides by the + # channel's stretch range, so they scale back up by the same range. + "cmy": tuple(cmy_offsets[ch] * density_range for ch in range(3)), + "exposure_on": tuple(exposure_offset != 0.0 or cmy_offsets[ch] != 0.0 for ch in range(3)), + "contrast": contrast, + "contrast_on": contrast != 1.0, + "pivot": float(c["transfer_contrast_pivot"]), + "shadow_cmy": tuple(shadow_cmy), + "highlight_cmy": tuple(highlight_cmy), + "regional_on": tuple(shadow_cmy[ch] != 0.0 or highlight_cmy[ch] != 0.0 for ch in range(3)), + "wb_c": wb_c, + "wb_k": wb_k, + "shadow_density": shadow_density, + "highlight_density": highlight_density, + "zone_on": shadow_density != 0.0 or highlight_density != 0.0, + "zone_sh": sh_c, + "zone_hi": hi_c, + "zone_k": zone_k, + "black_taper": ZONE_BLACK_TAPER, + "toe": tuple(toe), + "shoulder": tuple(shoulder), + "toe_width": tuple(toe_widths or (base_width,) * 3), + "shoulder_width": tuple(shoulder_widths or (base_width,) * 3), + "toe_knee": float(c["transfer_toe_knee"]), + "shoulder_knee": float(c["transfer_shoulder_knee"]), + "sep_k": sep_k3, + "sep_on": sep_k3 != (1.0, 1.0, 1.0), + "cap_level": SEPARATION_CAP_LEVEL, + "cap_offset": SEPARATION_CAP_SPREAD - SEPARATION_CAP_LEVEL, + "cap_softness": SEPARATION_CAP_SOFTNESS, + "damping": damping, + "damping_on": damping > 0.0, + "ref_spread": float(EXPOSURE_CONSTANTS["separation_damping_ref_spread"]), + "positive_source": bool(positive_source), + "gain": 2.0 ** float(c["transfer_baseline_ev"]), + } + return ensure_image(negpy_core.apply_transfer_curve(np.ascontiguousarray(img_norm, dtype=np.float32), params)) def transfer_point_offsets(process: ProcessConfig) -> Tuple[Tuple[float, float, float], Tuple[float, float, float]]: diff --git a/negpy/features/transparency/shaders/transfer.wgsl b/negpy/features/transparency/shaders/transfer.wgsl index 2047ac389..b360b96ab 100644 --- a/negpy/features/transparency/shaders/transfer.wgsl +++ b/negpy/features/transparency/shaders/transfer.wgsl @@ -40,7 +40,7 @@ struct TransferUniforms { @group(0) @binding(1) var output_tex: texture_storage_2d; @group(0) @binding(2) var params: TransferUniforms; -// width * log(1 + exp(x / width)), overflow-safe — mirrors logic.py::_softplus. +// width * log(1 + exp(x / width)), overflow-safe — mirrors softplus in negpy_core/src/transfer.rs. fn softplus(x: f32, width: f32) -> f32 { let t = x / width; return width * (log(1.0 + exp(-abs(t))) + max(t, 0.0)); diff --git a/negpy_core/src/lib.rs b/negpy_core/src/lib.rs index 6a7b175a2..1f0206c45 100644 --- a/negpy_core/src/lib.rs +++ b/negpy_core/src/lib.rs @@ -43,6 +43,7 @@ pub mod skew; pub mod stats; pub mod stitch; pub mod toning; +pub mod transfer; pub mod warp; pub mod zones; @@ -83,6 +84,7 @@ fn negpy_core(m: &pyo3::Bound<'_, pyo3::types::PyModule>) -> pyo3::PyResult<()> py::bounds::register(m)?; py::color::register(m)?; py::toning::register(m)?; + py::transfer::register(m)?; py::sensor::register(m)?; py::hdr::register(m)?; py::display::register(m)?; diff --git a/negpy_core/src/py/mod.rs b/negpy_core/src/py/mod.rs index e4caf415d..5c59c3ed4 100644 --- a/negpy_core/src/py/mod.rs +++ b/negpy_core/src/py/mod.rs @@ -33,6 +33,7 @@ pub mod skew; pub mod stats; pub mod stitch; pub mod toning; +pub mod transfer; pub mod warp; pub mod zones; diff --git a/negpy_core/src/py/transfer.rs b/negpy_core/src/py/transfer.rs new file mode 100644 index 000000000..8eff44cd1 --- /dev/null +++ b/negpy_core/src/py/transfer.rs @@ -0,0 +1,98 @@ +use numpy::{PyArrayDyn, PyArrayMethods, PyReadonlyArrayDyn, PyUntypedArrayMethods}; +use pyo3::exceptions::PyValueError; +use pyo3::prelude::*; +use pyo3::types::PyDict; + +use crate::transfer::{self, Transfer}; + +fn get<'py, T: FromPyObjectOwned<'py>>(p: &Bound<'py, PyDict>, key: &str) -> PyResult { + p.get_item(key)? + .ok_or_else(|| PyValueError::new_err(format!("missing transfer parameter {key}")))? + .extract() + .map_err(Into::into) +} + +fn f3<'py>(p: &Bound<'py, PyDict>, key: &str) -> PyResult<[f32; 3]> { + let v: [f64; 3] = get(p, key)?; + Ok(v.map(|x| x as f32)) +} + +fn b3<'py>(p: &Bound<'py, PyDict>, key: &str) -> PyResult<[bool; 3]> { + get(p, key) +} + +fn f1<'py>(p: &Bound<'py, PyDict>, key: &str) -> PyResult { + Ok(get::(p, key)? as f32) +} + +/// `apply_transfer_curve(img, params)` of a contiguous HxWx3 float32 image; `params` is the dict +/// `transparency.logic._core_params` builds, every value as NumPy would round it to float32. +#[pyfunction] +fn apply_transfer_curve<'py>( + py: Python<'py>, + img: PyReadonlyArrayDyn<'py, f32>, + params: &Bound<'py, PyDict>, +) -> PyResult>> { + let [_, _, 3] = img.shape()[..] else { + return Err(PyValueError::new_err("needs an HxWx3 image")); + }; + let p = params; + let t = Transfer { + density_range: f1(p, "density_range")?, + cast_gain: f3(p, "cast_gain")?, + cast_offset: f3(p, "cast_offset")?, + cast_on: b3(p, "cast_on")?, + exposure_offset: f1(p, "exposure_offset")?, + cmy: f3(p, "cmy")?, + exposure_on: b3(p, "exposure_on")?, + contrast: f1(p, "contrast")?, + contrast_on: get(p, "contrast_on")?, + pivot: f1(p, "pivot")?, + shadow_cmy: f3(p, "shadow_cmy")?, + highlight_cmy: f3(p, "highlight_cmy")?, + regional_on: b3(p, "regional_on")?, + wb_c: f1(p, "wb_c")?, + wb_k: f1(p, "wb_k")?, + shadow_density: f1(p, "shadow_density")?, + highlight_density: f1(p, "highlight_density")?, + zone_on: get(p, "zone_on")?, + zone_sh: f1(p, "zone_sh")?, + zone_hi: f1(p, "zone_hi")?, + zone_k: f1(p, "zone_k")?, + black_taper: f1(p, "black_taper")?, + toe: f3(p, "toe")?, + shoulder: f3(p, "shoulder")?, + toe_width: f3(p, "toe_width")?, + shoulder_width: f3(p, "shoulder_width")?, + toe_knee: f1(p, "toe_knee")?, + shoulder_knee: f1(p, "shoulder_knee")?, + sep_k: f3(p, "sep_k")?, + sep_on: get(p, "sep_on")?, + cap_level: f1(p, "cap_level")?, + cap_offset: f1(p, "cap_offset")?, + cap_softness: f1(p, "cap_softness")?, + damping: f1(p, "damping")?, + damping_on: get(p, "damping_on")?, + ref_spread: f1(p, "ref_spread")?, + positive_source: get(p, "positive_source")?, + gain: f1(p, "gain")?, + }; + let src = img + .as_slice() + .map_err(|_| PyValueError::new_err("needs a contiguous image"))?; + let out = PyArrayDyn::::zeros(py, img.shape(), false); + { + let mut rw = out.readwrite(); + let dst = rw + .as_slice_mut() + .map_err(|e| PyValueError::new_err(e.to_string()))?; + super::detach(py, dst.len(), || { + transfer::apply_transfer_curve(src, &t, dst) + }); + } + Ok(out) +} + +pub fn register(m: &Bound<'_, PyModule>) -> PyResult<()> { + m.add_function(wrap_pyfunction!(apply_transfer_curve, m)?) +} diff --git a/negpy_core/src/transfer.rs b/negpy_core/src/transfer.rs new file mode 100644 index 000000000..9acc68ca8 --- /dev/null +++ b/negpy_core/src/transfer.rs @@ -0,0 +1,240 @@ +//! `apply_transfer_curve` of negpy/features/transparency/logic.py: normalized log density to a +//! scene-linear positive, the render of every slide. Each step keeps NumPy's float32 order; the +//! transcendentals run in float64 (`mathf`) and round once, so they sit within a few ULP of +//! NumPy's SIMD float32 functions and give the same bits on every platform. + +use crate::mathf::{exp_f64, exp10_f64, log1p_f64, pow_f64, tanh_f64}; +use rayon::prelude::*; + +/// Every scalar the curve reads, already in float32 as NumPy's NEP 50 rule rounds them. +#[derive(Clone, Copy)] +pub struct Transfer { + pub density_range: f32, + pub cast_gain: [f32; 3], + pub cast_offset: [f32; 3], + pub cast_on: [bool; 3], + pub exposure_offset: f32, + pub cmy: [f32; 3], + pub exposure_on: [bool; 3], + pub contrast: f32, + pub contrast_on: bool, + pub pivot: f32, + pub shadow_cmy: [f32; 3], + pub highlight_cmy: [f32; 3], + pub regional_on: [bool; 3], + pub wb_c: f32, + pub wb_k: f32, + pub shadow_density: f32, + pub highlight_density: f32, + pub zone_on: bool, + pub zone_sh: f32, + pub zone_hi: f32, + pub zone_k: f32, + pub black_taper: f32, + pub toe: [f32; 3], + pub shoulder: [f32; 3], + pub toe_width: [f32; 3], + pub shoulder_width: [f32; 3], + pub toe_knee: f32, + pub shoulder_knee: f32, + pub sep_k: [f32; 3], + pub sep_on: bool, + pub cap_level: f32, + pub cap_offset: f32, + pub cap_softness: f32, + pub damping: f32, + pub damping_on: bool, + pub ref_spread: f32, + pub positive_source: bool, + pub gain: f32, +} + +/// `np.maximum(a, b)` with b not NaN: NaN in `a` propagates. +#[inline(always)] +fn np_max(a: f32, b: f32) -> f32 { + if a.is_nan() || a >= b { a } else { b } +} + +/// `_sigmoid`: 0.5 * (1 + tanh(0.5 x)) in float32. +#[inline(always)] +fn sigmoid(x: f32) -> f32 { + 0.5f32 * (1.0f32 + tanh_f64((0.5f32 * x) as f64) as f32) +} + +/// `_softplus`: width * (logaddexp(0, -|t|) + max(t, 0)) with t = x / width. +#[inline(always)] +fn softplus(x: f32, width: f32) -> f32 { + let t = x / width; + let lae = log1p_f64(exp_f64(-(t.abs() as f64))) as f32; + width * (lae + np_max(t, 0.0)) +} + +/// Pixels per staged run: each step is a flat loop over a run, so the float64 functions vectorize. +const RUN: usize = 256; + +impl Transfer { + /// `shape(d, ch)` over a run: one channel's capture density through every control above + /// Dye Separation. + #[inline(always)] + fn shape(&self, d: &mut [f32], ch: usize) { + if self.cast_on[ch] { + let (g, o) = (self.cast_gain[ch], self.cast_offset[ch]); + d.iter_mut().for_each(|v| *v = *v * g + o); + } + if self.exposure_on[ch] { + let (e, c) = (self.exposure_offset, self.cmy[ch]); + d.iter_mut().for_each(|v| *v = *v - e + c); + } + if self.contrast_on { + let (p, k) = (self.pivot, self.contrast); + d.iter_mut().for_each(|v| *v = p + (*v - p) * k); + } + if self.regional_on[ch] { + let (k, c, s, h) = ( + self.wb_k, + self.wb_c, + self.shadow_cmy[ch], + self.highlight_cmy[ch], + ); + d.iter_mut().for_each(|v| { + let w_sh = sigmoid(k * (*v - c)); + *v = *v + s * w_sh + h * (1.0f32 - w_sh); + }); + } + if self.zone_on { + let t = *self; + d.iter_mut().for_each(|v| { + let x = *v; + let b = ((t.density_range - x) / t.black_taper).clamp(0.0, 1.0); + let w_sh = sigmoid(t.zone_k * (x - t.zone_sh)) * (b * b * (3.0f32 - 2.0f32 * b)); + let w_hi = 1.0f32 - sigmoid(t.zone_k * (x - t.zone_hi)); + *v = x + t.shadow_density * w_sh + t.highlight_density * w_hi; + }); + } + if self.toe[ch] != 0.0 { + let (a, k, w) = (self.toe[ch], self.toe_knee, self.toe_width[ch]); + d.iter_mut().for_each(|v| *v -= a * softplus(*v - k, w)); + } + if self.shoulder[ch] != 0.0 { + let (a, k, w) = ( + self.shoulder[ch], + self.shoulder_knee, + self.shoulder_width[ch], + ); + d.iter_mut().for_each(|v| *v += a * softplus(k - *v, w)); + } + } + + /// Dye Separation over a run: `dens` (the shaped channels) pulled from or toward the + /// reference, the mean of the capped and shaped channels. + #[inline(always)] + fn separate(&self, capture: &[[f32; RUN]; 3], dens: &mut [[f32; RUN]; 3], n: usize) { + let (soft, level, off) = (self.cap_softness, self.cap_level, self.cap_offset); + let mut cap = [0f32; RUN]; + for i in 0..n { + let (a, b, c) = (capture[0][i], capture[1][i], capture[2][i]); + let lo = if a.is_nan() || b.is_nan() || c.is_nan() { + f32::NAN + } else { + a.min(b).min(c) + }; + cap[i] = level + softplus(lo + off, soft); + } + let mut capped = [[0f32; RUN]; 3]; + for ch in 0..3 { + for i in 0..n { + capped[ch][i] = cap[i] - softplus(cap[i] - capture[ch][i], soft); + } + self.shape(&mut capped[ch][..n], ch); + } + for i in 0..n { + let d_ref = (capped[0][i] + capped[1][i] + capped[2][i]) / 3.0f32; + let e: [f32; 3] = std::array::from_fn(|ch| dens[ch][i] - d_ref); + if self.damping_on { + let ce: [f32; 3] = std::array::from_fn(|ch| capped[ch][i] - d_ref); + let (a, b, c) = (ce[0] - ce[1], ce[1] - ce[2], ce[0] - ce[2]); + let chroma = ((a * a + b * b + c * c) / 3.0f32).sqrt(); + let h = (self.ref_spread - chroma) / (self.ref_spread + chroma); + let expo = (1.0f32 - self.damping) + self.damping * h; + for ch in 0..3 { + let k = self.sep_k[ch]; + let kf = pow_f64(k as f64, expo as f64) as f32; + let kf = if kf.is_nan() || kf <= 3.0 { kf } else { 3.0 }; + dens[ch][i] = d_ref + if k <= 0.0 { 0.0 } else { kf } * e[ch]; + } + } else { + for ch in 0..3 { + dens[ch][i] = d_ref + self.sep_k[ch] * e[ch]; + } + } + } + } +} + +#[multiversion::multiversion(targets( + "x86_64+avx512f+avx512bw+avx512vl+avx512dq", + "x86_64+avx2+fma", + "aarch64+neon" +))] +fn transfer_block(t: &Transfer, src: &[[f32; 3]], dst: &mut [[f32; 3]]) { + let mut capture = [[0f32; RUN]; 3]; + let mut dens = [[0f32; RUN]; 3]; + for (o, p) in dst.chunks_mut(RUN).zip(src.chunks(RUN)) { + let n = p.len(); + for ch in 0..3 { + for i in 0..n { + capture[ch][i] = p[i][ch] * t.density_range; + } + dens[ch][..n].copy_from_slice(&capture[ch][..n]); + t.shape(&mut dens[ch][..n], ch); + } + if t.sep_on { + t.separate(&capture, &mut dens, n); + } + for ch in 0..3 { + for i in 0..n { + let v = exp10_f64(-(dens[ch][i] as f64)) as f32; + o[i][ch] = if t.positive_source { + v.clamp(0.0, 1.0) + } else { + let x = np_max(v * t.gain, 0.0); + let num = x * (2.51f32 * x + 0.03f32); + let den = x * (2.43f32 * x + 0.59f32) + 0.14f32; + (num / np_max(den, 1e-8)).clamp(0.0, 1.0) + }; + } + } + } +} + +/// The transfer curve of an HxWx3 normalized image into `dst`. +pub fn apply_transfer_curve(src: &[f32], t: &Transfer, dst: &mut [f32]) { + let t = *t; + dst.par_chunks_mut(3 * 4096) + .zip(src.par_chunks(3 * 4096)) + .for_each(move |(d, s)| transfer_block(&t, s.as_chunks::<3>().0, d.as_chunks_mut::<3>().0)); +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn softplus_is_relu_far_from_the_knee() { + assert!((softplus(10.0, 0.45) - 10.0).abs() < 1e-5); + assert!(softplus(-10.0, 0.45) < 1e-6); + assert!((softplus(0.0, 1.0) - std::f32::consts::LN_2).abs() < 1e-6); + } + + #[test] + fn sigmoid_is_half_at_zero() { + assert_eq!(sigmoid(0.0), 0.5); + assert!(sigmoid(40.0) > 0.999_99); + } + + #[test] + fn np_max_keeps_nan() { + assert!(np_max(f32::NAN, 0.0).is_nan()); + assert_eq!(np_max(-1.0, 0.0), 0.0); + } +} diff --git a/tests/test_core_same_bits.py b/tests/test_core_same_bits.py index 4176a62e1..8f6381a60 100644 --- a/tests/test_core_same_bits.py +++ b/tests/test_core_same_bits.py @@ -6,6 +6,7 @@ from negpy.features.exposure.logic import apply_characteristic_curve from negpy.features.toning.logic import apply_chemical_toning +from negpy.features.transparency.logic import apply_transfer_curve from negpy.kernel.image.logic import WORKING_LAB # Recorded on Linux x86-64, where the cv2 oracles pass. Every release platform must give these @@ -75,6 +76,19 @@ def _filled(fill): "matrix_trc_u8": lambda: negpy_core.matrix_trc_u8((_IMG3 * 255).astype(np.uint8), np.ascontiguousarray(_RAMP[:, ::256]), _M, _RAMP), "vignette": lambda: negpy_core.apply_vignette(_IMG3, 0.7, 0.6, 0.3), "flat_curve": lambda: negpy_core.flat_curve(_IMG3, 0.65, 0.1, (0.013, -0.02, 0.0), _M), + "transfer_curve": lambda: apply_transfer_curve( + _IMG3, + 0.1, + 1.1, + (0.2, 0.3, 0.1), + (0.2, 0.0, 0.25), + (0.01, 0.0, -0.02), + shadow_density=0.15, + shadow_cmy=(0.02, 0.0, 0.0), + cast_gain=(1.02, 1.0, 0.97), + separation=1.3, + damping=0.5, + ), "matrix_trc_u16": lambda: negpy_core.matrix_trc_u16( (_IMG3 * 65535).astype(np.uint16), np.ascontiguousarray(np.tile(_LIN, (3, 1))), _M, np.ascontiguousarray(_RAMP[:, ::16]) ), @@ -148,6 +162,7 @@ def _filled(fill): "matrix_trc_u16": "8fd1af40718fc9bc59541e1128e46cb32e08ea05626f7cbd4cebe7905cecd02a", "vignette": "dea377797d6b0852a8dec1665c95aa87c9ab2b6a8e03f808e27140bc4b6566c6", "flat_curve": "e47f6084a58e19119c6e3364ebf017dfc7f039ec6c73e0eea8e942519021bf0b", + "transfer_curve": "afdb5891162dca5abe4530bd098469b26a0be16a41df32f4f2928ddb8d37b943", "usm_sharpen": "7267d2a3976c9381eadde01fe11ea60d46345b4747283e6ac68a2aa4a025d0a8", "sabattier": "f8956f67786c7c1d0e628af4b82fb12b6dde42e70da2504585f5bf20e7fbb461", "lith": "98dc05c62e191da2ad00c1f2327c48beff5af7e882be64cac5302cb0ea5099e5", From 1c9e5345af9ef0dce1cea4859e7e8b31046bd861 Mon Sep 17 00:00:00 2001 From: Marcin Zawalski Date: Sun, 25 Oct 2026 09:51:38 +0100 Subject: [PATCH 09/12] perf(loaders): TIFF and DNG preview strips in Rust A bounded TIFF preview mapped every decoded strip through the display LUT in NumPy (a full-resolution gather) and shrank it with PIL's BOX resize; the LinearRaw DNG preview ran its whole display curve, np.power included, on every tile. negpy_core.box_into_u8 area-averages a strip straight into its preview rows, reading each sample through the display table. The DNG curve becomes one 3 x 65536 table per file, built with the same formula. Previews stay within 2 levels of PIL's BOX output. 23.5 MP deflate TIFF to a 512 px preview: 338 -> 226 ms per frame; what is left is strip inflate on the decode workers. --- negpy/infrastructure/loaders/helpers.py | 50 +++++++++-------- negpy_core/src/py/resize.rs | 75 +++++++++++++++++++++++++ negpy_core/src/resize.rs | 72 ++++++++++++++++++++++++ 3 files changed, 174 insertions(+), 23 deletions(-) diff --git a/negpy/infrastructure/loaders/helpers.py b/negpy/infrastructure/loaders/helpers.py index f6325bdcc..404b83d34 100644 --- a/negpy/infrastructure/loaders/helpers.py +++ b/negpy/infrastructure/loaders/helpers.py @@ -4,6 +4,7 @@ from types import SimpleNamespace from typing import Any, Callable, Optional, Tuple +import negpy_core import numpy as np import rawpy from PIL import Image, ImageCms @@ -251,6 +252,7 @@ def _reduced_top_level_page(tif: Any) -> Optional[Any]: _linear_u16_levels[_linear_u16_low] *= 4.5 _linear_u16_levels[~_linear_u16_low] = 1.099 * np.power(_linear_u16_levels[~_linear_u16_low], 1.0 / 2.222) - 0.099 _LINEAR_U16_DISPLAY_LUT = np.clip(_linear_u16_levels * 255.0, 0, 255).astype(np.uint8) +_LINEAR_U16_DISPLAY_ROWS = np.ascontiguousarray(_LINEAR_U16_DISPLAY_LUT[None, :]) del _linear_u16_levels, _linear_u16_low @@ -362,6 +364,26 @@ def dng_quick_preview(file_path: str) -> Optional[Image.Image]: return Image.fromarray(np.ascontiguousarray(arr)) +def _dng_display_lut(dtype: Any, linearization: Optional[np.ndarray], black: np.ndarray, white: np.ndarray, wb: np.ndarray) -> np.ndarray: + """The LinearRaw preview's display value of every code, per channel (3 x 65536 uint8): + linearization, black/white normalize, as-shot white balance, the BT.709 curve.""" + codes = np.arange(int(np.iinfo(dtype).max) + 1, dtype=dtype) + source = np.repeat(codes[:, None, None], 3, axis=2) + if linearization is not None: + indices = np.clip(source, 0, linearization.size - 1).astype(np.int32) + data = linearization[indices] + else: + data = source.astype(np.float32) + data = np.clip((data - black) / np.maximum(white - black, 1e-6), 0.0, 1.0) + data = np.clip(data * wb, 0.0, 1.0) + low = data < 0.018 + data[low] *= 4.5 + data[~low] = 1.099 * np.power(data[~low], 1.0 / 2.222) - 0.099 + table = np.zeros((3, 65536), dtype=np.uint8) + table[:, : codes.size] = np.clip(data * 255.0, 0, 255).astype(np.uint8)[:, 0, :].T + return table + + def dng_bounded_preview( file_path: str, max_edge: int, @@ -424,6 +446,7 @@ def tag(name: str) -> Any: crop_size = _dng_tag_floats(tag("DefaultCropSize")) orientation_tag = tag("Orientation") orientation = int(orientation_tag.value) if orientation_tag is not None else 1 + lut = _dng_display_lut(page.dtype, linearization, black, white, wb) for decoded, position, _shape in page.segments(maxworkers=1): if should_cancel is not None and should_cancel(): raise InterruptedError("thumbnail cancelled") @@ -437,27 +460,14 @@ def tag(name: str) -> Any: valid_width = min(tile.shape[1], width - x) if valid_height <= 0 or valid_width <= 0: continue - source = tile[:valid_height, :valid_width, :3] - if linearization is not None: - indices = np.clip(source, 0, linearization.size - 1).astype(np.int32) - data = linearization[indices] - else: - data = source.astype(np.float32) - data = np.clip((data - black) / np.maximum(white - black, 1e-6), 0.0, 1.0) - data = np.clip(data * wb, 0.0, 1.0) - low = data < 0.018 - data[low] *= 4.5 - data[~low] = 1.099 * np.power(data[~low], 1.0 / 2.222) - 0.099 - tile_u8 = np.clip(data * 255.0, 0, 255).astype(np.uint8) - left = int(round(x * out_width / width)) top = int(round(y * out_height / height)) right = int(round((x + valid_width) * out_width / width)) bottom = int(round((y + valid_height) * out_height / height)) if right <= left or bottom <= top: continue - small = Image.fromarray(tile_u8).resize((right - left, bottom - top), Image.Resampling.BOX) - output[top:bottom, left:right] = np.asarray(small) + source = np.ascontiguousarray(tile[:valid_height, :valid_width, :3], dtype=np.uint16) + negpy_core.box_into_u8(source, output[top:bottom, left:right], lut) image = Image.fromarray(output) if crop_origin.size >= 2 and crop_size.size >= 2: @@ -521,20 +531,14 @@ def bounded_tiff_page_preview( output = np.zeros((out_height, out_width, 3), dtype=np.uint8) def place(source: np.ndarray, y: int, x: int) -> None: - if source.dtype == np.uint16: - source = linear_uint16_to_display_uint8(source) - if source.shape[2] == 1: - source = np.repeat(source, 3, axis=2) - elif source.shape[2] == 4: - source = source[:, :, :3] left = int(round(x * out_width / width)) top = int(round(y * out_height / height)) right = int(round((x + source.shape[1]) * out_width / width)) bottom = int(round((y + source.shape[0]) * out_height / height)) if right <= left or bottom <= top: return - small = Image.fromarray(np.ascontiguousarray(source)).resize((right - left, bottom - top), Image.Resampling.BOX) - output[top:bottom, left:right] = np.asarray(small) + lut = _LINEAR_U16_DISPLAY_ROWS if source.dtype == np.uint16 else None + negpy_core.box_into_u8(np.ascontiguousarray(source), output[top:bottom, left:right], lut) if banded: # An uncompressed page is one run of rows at its first offset, so it reads in diff --git a/negpy_core/src/py/resize.rs b/negpy_core/src/py/resize.rs index 6eb46078b..adf6d1226 100644 --- a/negpy_core/src/py/resize.rs +++ b/negpy_core/src/py/resize.rs @@ -336,7 +336,82 @@ fn area_rows_into( Ok(()) } +fn box_into>( + py: Python<'_>, + src: &PyReadonlyArrayDyn<'_, T>, + dst: &mut numpy::PyReadwriteArray3<'_, u8>, + lut: Option<(&[u8], usize)>, +) -> PyResult<()> { + let (h, w, cn) = match src.shape()[..] { + [h, w] => (h, w, 1), + [h, w, c] if c >= 1 => (h, w, c), + _ => { + return Err(PyValueError::new_err( + "box_into needs an HxW or HxWxC source", + )); + } + }; + let s = src + .as_slice() + .map_err(|_| PyValueError::new_err("needs a contiguous source"))?; + let (oh, ow, oc) = dst.as_array().dim(); + if oc != 3 || h == 0 || w == 0 { + return Err(PyValueError::new_err("box_into needs an HxWx3 destination")); + } + let mut buf = vec![0u8; oh * ow * 3]; + super::detach(py, 0, || { + resize::box_into_u8(s, h, w, cn, lut, &mut buf, oh, ow, ow * 3) + }); + for (mut row, packed) in dst + .as_array_mut() + .outer_iter_mut() + .zip(buf.chunks_exact(ow * 3)) + { + for (mut px, p) in row.outer_iter_mut().zip(packed.as_chunks::<3>().0) { + for c in 0..3 { + px[c] = p[c]; + } + } + } + Ok(()) +} + +/// `box_into(src, dst, lut=None)`: a uint8 or uint16 image area-averaged into the uint8 HxWx3 +/// view `dst`, uint16 samples read through `lut` (k x 65536 uint8, one row per channel or one). +#[pyfunction] +#[pyo3(signature = (src, dst, lut=None))] +fn box_into_u8( + py: Python<'_>, + src: &Bound<'_, PyAny>, + mut dst: numpy::PyReadwriteArray3<'_, u8>, + lut: Option>, +) -> PyResult<()> { + let table = match &lut { + Some(t) => { + let k = t.shape()[0]; + if t.shape()[1] != 65536 || k == 0 { + return Err(PyValueError::new_err("the LUT must be k x 65536")); + } + Some(( + t.as_slice() + .map_err(|_| PyValueError::new_err("needs a contiguous LUT"))?, + k, + )) + } + None => None, + }; + if let Ok(a) = src.extract::>() { + if table.is_none() { + return Err(PyValueError::new_err("a uint16 source needs a LUT")); + } + return box_into(py, &a, &mut dst, table); + } + let a = src.extract::>()?; + box_into(py, &a, &mut dst, None) +} + pub fn register(m: &Bound<'_, PyModule>) -> PyResult<()> { + m.add_function(wrap_pyfunction!(box_into_u8, m)?)?; m.add_function(wrap_pyfunction!(area_rows_into, m)?)?; m.add_function(wrap_pyfunction!(resize_rows_linear, m)?)?; m.add_function(wrap_pyfunction!(resize_u8, m)?)?; diff --git a/negpy_core/src/resize.rs b/negpy_core/src/resize.rs index ea8a4ffac..ed829e482 100644 --- a/negpy_core/src/resize.rs +++ b/negpy_core/src/resize.rs @@ -1051,10 +1051,82 @@ pub fn area_rows( }); } +/// The overlap of source pixels with output pixel `o` of `n_out` along an axis of `n_in`: +/// (first source index, weights summing to 1). +fn box_weights(o: usize, n_in: usize, n_out: usize) -> (usize, Vec) { + let scale = n_in as f64 / n_out as f64; + let (a, b) = (o as f64 * scale, (o + 1) as f64 * scale); + let first = a.floor() as usize; + let last = (b.ceil() as usize).min(n_in); + let mut w: Vec = (first..last) + .map(|i| ((i + 1) as f64).min(b) - (i as f64).max(a)) + .collect(); + let sum: f64 = w.iter().sum(); + w.iter_mut().for_each(|v| *v /= sum); + (first, w) +} + +/// An (h, w, cn) uint8 or uint16 image area-averaged into the (oh, ow, 3) uint8 `dst`, each +/// sample read through `lut` (`k` tables of 65536 codes, one per channel or one shared) when +/// given. One channel fills all three; a fourth is dropped. Rounded to nearest. +#[allow(clippy::too_many_arguments)] +pub fn box_into_u8>( + src: &[T], + h: usize, + w: usize, + cn: usize, + lut: Option<(&[u8], usize)>, + dst: &mut [u8], + oh: usize, + ow: usize, + dst_row: usize, +) { + let cols: Vec<(usize, Vec)> = (0..ow).map(|x| box_weights(x, w, ow)).collect(); + let value = |y: usize, x: usize, c: usize| -> f64 { + let sc = c.min(cn - 1); + let v: u32 = src[(y * w + x) * cn + sc].into(); + match lut { + Some((t, k)) => t[(c.min(k - 1)) * 65536 + v as usize] as f64, + None => v as f64, + } + }; + for oy in 0..oh { + let (y0, wy) = box_weights(oy, h, oh); + let row = &mut dst[oy * dst_row..oy * dst_row + ow * 3]; + for (ox, (x0, wx)) in cols.iter().enumerate() { + let mut acc = [0f64; 3]; + for (dy, &fy) in wy.iter().enumerate() { + for (dx, &fx) in wx.iter().enumerate() { + let f = fy * fx; + for (c, a) in acc.iter_mut().enumerate() { + *a += f * value(y0 + dy, x0 + dx, c); + } + } + } + for c in 0..3 { + row[ox * 3 + c] = acc[c].round().clamp(0.0, 255.0) as u8; + } + } + } +} + #[cfg(test)] mod tests { use super::*; + #[test] + fn box_into_u8_averages_through_the_lut() { + let src: Vec = vec![0, 2, 4, 6]; + let lut: Vec = (0..65536).map(|v| (v * 10).min(255) as u8).collect(); + let mut dst = [0u8; 3]; + box_into_u8(&src, 2, 2, 1, Some((&lut, 1)), &mut dst, 1, 1, 3); + assert_eq!(dst, [30, 30, 30]); + let rgb: Vec = vec![10, 20, 30, 255, 30, 40, 50, 255]; + let mut out = [0u8; 3]; + box_into_u8(&rgb, 1, 2, 4, None, &mut out, 1, 1, 3); + assert_eq!(out, [20, 30, 40]); + } + #[test] fn flat_image_stays_flat() { let src = vec![0.25f32; 37 * 53 * 3]; From 045654de82c97fd1a286d4864632e4498f352530 Mon Sep 17 00:00:00 2001 From: Marcin Zawalski Date: Sun, 25 Oct 2026 09:58:24 +0100 Subject: [PATCH 10/12] perf(export): parallel PNG encoder in Rust PNG exports went through PIL (8-bit) and imagecodecs/libpng (16-bit), both filtering and deflating the whole image on one thread. negpy_core.encode_png filters bands of 64 rows (the minimum-sum-of-absolute-differences choice over the five PNG filters) and deflates each band with miniz_oxide on its own core; sync flushes join the bands into one zlib stream and the Adler-32 is combined from theirs. encode_png splices the iCCP, pHYs, eXIf and iTXt chunks for both depths. Pixels decode identically in libpng and PIL. New crates: miniz_oxide, adler2, crc32fast (pure Rust; NOTICE.md). 23.5 MP, level 6, 24 threads: 8-bit 16.4 s -> 0.74 s (size +0.8%), 16-bit 16.4 s -> 1.1 s (size +0.05%). --- .github/workflows/ci.yml | 2 +- NOTICE.md | 3 +- negpy/services/export/encoders.py | 69 +++-------- negpy_core/Cargo.lock | 33 +++++ negpy_core/Cargo.toml | 3 + negpy_core/src/lib.rs | 1 + negpy_core/src/png.rs | 198 ++++++++++++++++++++++++++++++ negpy_core/src/py/jpeg.rs | 47 +++++++ tests/test_core_png.py | 35 ++++++ 9 files changed, 339 insertions(+), 52 deletions(-) create mode 100644 negpy_core/src/png.rs create mode 100644 tests/test_core_png.py diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 600acce41..06713880b 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -129,4 +129,4 @@ jobs: run: uv sync --group dev --group plustek - name: negpy_core pins and oracles - run: uv run pytest tests/test_percentiles.py tests/test_core_same_bits.py tests/test_core_orient.py tests/test_core_ir_downsample.py tests/test_core_flat_curve.py tests/test_core_resize.py tests/test_core_filter.py tests/test_core_morph.py tests/test_core_gray.py tests/test_core_contours.py tests/test_core_shapes.py tests/test_core_components.py tests/test_core_color.py tests/test_core_warp.py tests/test_core_lab.py tests/test_core_altprocess.py + run: uv run pytest tests/test_percentiles.py tests/test_core_same_bits.py tests/test_core_orient.py tests/test_core_ir_downsample.py tests/test_core_flat_curve.py tests/test_core_png.py tests/test_core_resize.py tests/test_core_filter.py tests/test_core_morph.py tests/test_core_gray.py tests/test_core_contours.py tests/test_core_shapes.py tests/test_core_components.py tests/test_core_color.py tests/test_core_warp.py tests/test_core_lab.py tests/test_core_altprocess.py diff --git a/NOTICE.md b/NOTICE.md index d16794516..0876ce521 100644 --- a/NOTICE.md +++ b/NOTICE.md @@ -65,7 +65,8 @@ Apache-2.0), ndarray, num-traits, num-complex, num-integer, matrixmultiply and rawpointer (MIT OR Apache-2.0), libc, once_cell, either, portable-atomic and rustc-hash (MIT OR Apache-2.0), multiversion and target-features (MIT OR Apache-2.0), libm (MIT, a port of musl's math library), rustfft, realfft, primal-check, strength_reduce and transpose (MIT OR -Apache-2.0). +Apache-2.0), miniz_oxide (MIT OR Zlib OR Apache-2.0), adler2 (0BSD OR MIT OR Apache-2.0), +crc32fast and cfg-if (MIT OR Apache-2.0). negpy_core/Cargo.lock lists the exact versions. negpy_core/src/resize.rs and negpy_core/src/filter.rs port OpenCV's float32 area diff --git a/negpy/services/export/encoders.py b/negpy/services/export/encoders.py index 73a28197e..35e8deabb 100644 --- a/negpy/services/export/encoders.py +++ b/negpy/services/export/encoders.py @@ -15,7 +15,7 @@ import negpy_core import numpy as np import tifffile -from PIL import Image, ImageFile, PngImagePlugin +from PIL import Image, ImageFile from negpy.domain.models import TiffCompression from negpy.features.metadata.resolution import Resolution @@ -67,28 +67,27 @@ def encode_png( exif: Optional[bytes] = None, xmp: Optional[bytes] = None, ) -> bytes: - """Encode uint8/uint16 pixels as PNG. - - uint16 goes through imagecodecs and gets its ancillary chunks spliced in, - because PIL has no 16-bit RGB mode. - """ - if arr.dtype == np.uint16 and arr.ndim == 3: - return _encode_png_u16(arr, icc=icc, resolution=resolution, level=level, exif=exif, xmp=xmp) - - buf = io.BytesIO() - kwargs: Dict[str, Any] = {"format": "PNG", "compress_level": level} - if resolution is not None: - kwargs["dpi"] = (resolution.x_dpi, resolution.y_dpi) + """Encode uint8/uint16 pixels as PNG: negpy_core filters and deflates bands of rows on + every core, and the ancillary chunks are spliced in ahead of the image data.""" + png = bytes(negpy_core.encode_png(np.ascontiguousarray(arr), int(level))) + chunks = [] if icc: - kwargs["icc_profile"] = icc + chunks.append(_png_chunk(b"iCCP", b"ICC Profile\x00\x00" + zlib.compress(icc, 9))) + if resolution is not None: + ppm_x = max(1, int(round(resolution.x_dpi * _INCH_PER_METRE))) + ppm_y = max(1, int(round(resolution.y_dpi * _INCH_PER_METRE))) + chunks.append(_png_chunk(b"pHYs", struct.pack(">IIB", ppm_x, ppm_y, _PNG_METRE))) if exif: - kwargs["exif"] = exif + # eXIf carries the bare TIFF structure; the APP1 marker piexif emits is a + # JPEG framing detail and makes readers reject the chunk. + chunks.append(_png_chunk(b"eXIf", exif[6:] if exif.startswith(b"Exif\x00\x00") else exif)) if xmp: - info = PngImagePlugin.PngInfo() - info.add_itxt("XML:com.adobe.xmp", xmp.decode("utf-8"), zip=False) - kwargs["pnginfo"] = info - Image.fromarray(arr).save(buf, **kwargs) - return buf.getvalue() + chunks.append(_png_chunk(b"iTXt", b"XML:com.adobe.xmp\x00\x00\x00\x00\x00" + xmp)) + if not chunks: + return png + # Ancillary chunks must precede the image data. + cut = png.index(b"IDAT") - 4 + return png[:cut] + b"".join(chunks) + png[cut:] def encode_jpeg( @@ -178,33 +177,3 @@ def encode_jxl( def _png_chunk(tag: bytes, data: bytes) -> bytes: return struct.pack(">I", len(data)) + tag + data + struct.pack(">I", zlib.crc32(tag + data) & 0xFFFFFFFF) - - -def _encode_png_u16( - arr: np.ndarray, - *, - icc: Optional[bytes], - resolution: Optional[Resolution], - level: int, - exif: Optional[bytes], - xmp: Optional[bytes], -) -> bytes: - png = bytes(imagecodecs.png_encode(np.ascontiguousarray(arr), level=level)) - chunks = [] - if icc: - chunks.append(_png_chunk(b"iCCP", b"ICC Profile\x00\x00" + zlib.compress(icc, 9))) - if resolution is not None: - ppm_x = max(1, int(round(resolution.x_dpi * _INCH_PER_METRE))) - ppm_y = max(1, int(round(resolution.y_dpi * _INCH_PER_METRE))) - chunks.append(_png_chunk(b"pHYs", struct.pack(">IIB", ppm_x, ppm_y, _PNG_METRE))) - if exif: - # eXIf carries the bare TIFF structure; the APP1 marker piexif emits is a - # JPEG framing detail and makes readers reject the chunk. - chunks.append(_png_chunk(b"eXIf", exif[6:] if exif.startswith(b"Exif\x00\x00") else exif)) - if xmp: - chunks.append(_png_chunk(b"iTXt", b"XML:com.adobe.xmp\x00\x00\x00\x00\x00" + xmp)) - if not chunks: - return png - # Ancillary chunks must precede the image data. - cut = png.index(b"IDAT") - 4 - return png[:cut] + b"".join(chunks) + png[cut:] diff --git a/negpy_core/Cargo.lock b/negpy_core/Cargo.lock index f13707c69..47f96063e 100644 --- a/negpy_core/Cargo.lock +++ b/negpy_core/Cargo.lock @@ -2,12 +2,33 @@ # 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" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "f2032f911046de80f0a198e0901378627c33f59ea0ac00e363d481118bd70a53" +[[package]] +name = "cfg-if" +version = "1.0.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4e7648175b45a9a48536d676f68d918270699102aa8dab5496df06904c914600" + +[[package]] +name = "crc32fast" +version = "1.5.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0ba106e844045ddc02c20768ac4b0b7cb45f2cc2ecbd833ac9f5dacb183c2055" +dependencies = [ + "cfg-if", +] + [[package]] name = "crossbeam-deque" version = "0.8.8" @@ -67,6 +88,15 @@ dependencies = [ "rawpointer", ] +[[package]] +name = "miniz_oxide" +version = "0.8.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1fa76a2c86f704bdb222d66965fb3d63269ce38518b83cb0575fca855ebb6316" +dependencies = [ + "adler2", +] + [[package]] name = "multiversion" version = "0.8.1" @@ -108,7 +138,10 @@ dependencies = [ name = "negpy_core" version = "0.1.0" dependencies = [ + "adler2", + "crc32fast", "libm", + "miniz_oxide", "multiversion", "numpy", "pyo3", diff --git a/negpy_core/Cargo.toml b/negpy_core/Cargo.toml index 70ab2811e..737a2d553 100644 --- a/negpy_core/Cargo.toml +++ b/negpy_core/Cargo.toml @@ -20,3 +20,6 @@ multiversion = "0.8" libm = "0.2" rustfft = "6" realfft = "3" +miniz_oxide = "0.8" +crc32fast = "1" +adler2 = "2" diff --git a/negpy_core/src/lib.rs b/negpy_core/src/lib.rs index 1f0206c45..3944b3391 100644 --- a/negpy_core/src/lib.rs +++ b/negpy_core/src/lib.rs @@ -29,6 +29,7 @@ pub mod meter; pub mod morph; pub mod orient; pub mod phase; +pub mod png; pub mod print_curve; pub mod process; pub mod register; diff --git a/negpy_core/src/png.rs b/negpy_core/src/png.rs new file mode 100644 index 000000000..0a9812200 --- /dev/null +++ b/negpy_core/src/png.rs @@ -0,0 +1,198 @@ +//! A PNG encoder that filters and deflates bands of rows in parallel. Each band is an +//! independent deflate run ended by a sync flush, so the bands join into one zlib stream; the +//! Adler-32 of the whole image is combined from the bands' own. + +use miniz_oxide::deflate::core::{ + CompressorOxide, TDEFLFlush, TDEFLStatus, compress_to_output, create_comp_flags_from_zip_params, +}; +use rayon::prelude::*; + +/// Rows per band. Each band restarts the deflate window, which costs a few bytes per band. +const BAND_ROWS: usize = 64; +const ADLER_BASE: u32 = 65521; + +/// zlib's adler32_combine: the checksum of A then B from theirs and B's length. +fn adler32_combine(a: u32, b: u32, len_b: usize) -> u32 { + let rem = (len_b % ADLER_BASE as usize) as u32; + let mut sum1 = a & 0xffff; + let mut sum2 = (rem * sum1) % ADLER_BASE; + sum1 += (b & 0xffff) + ADLER_BASE - 1; + sum2 += ((a >> 16) & 0xffff) + ((b >> 16) & 0xffff) + ADLER_BASE - rem; + if sum1 >= ADLER_BASE { + sum1 -= ADLER_BASE; + } + if sum1 >= ADLER_BASE { + sum1 -= ADLER_BASE; + } + if sum2 >= ADLER_BASE << 1 { + sum2 -= ADLER_BASE << 1; + } + if sum2 >= ADLER_BASE { + sum2 -= ADLER_BASE; + } + sum1 | (sum2 << 16) +} + +#[inline(always)] +fn paeth(a: u8, b: u8, c: u8) -> u8 { + let (ia, ib, ic) = (a as i16, b as i16, c as i16); + let p = ia + ib - ic; + let (pa, pb, pc) = ((p - ia).abs(), (p - ib).abs(), (p - ic).abs()); + if pa <= pb && pa <= pc { + a + } else if pb <= pc { + b + } else { + c + } +} + +/// One row with its filter byte, by the minimum sum of absolute differences over the five +/// PNG filters (the heuristic libpng suggests). +fn filter_row( + row: &[u8], + prev: Option<&[u8]>, + bpp: usize, + out: &mut Vec, + scratch: &mut [Vec; 5], +) { + let n = row.len(); + for s in scratch.iter_mut() { + s.clear(); + } + for i in 0..n { + let x = row[i]; + let a = if i >= bpp { row[i - bpp] } else { 0 }; + let b = prev.map_or(0, |p| p[i]); + let c = if i >= bpp { + prev.map_or(0, |p| p[i - bpp]) + } else { + 0 + }; + scratch[0].push(x); + scratch[1].push(x.wrapping_sub(a)); + scratch[2].push(x.wrapping_sub(b)); + scratch[3].push(x.wrapping_sub(((a as u16 + b as u16) / 2) as u8)); + scratch[4].push(x.wrapping_sub(paeth(a, b, c))); + } + let cost = |s: &Vec| { + s.iter() + .map(|&v| (v as i8).unsigned_abs() as u64) + .sum::() + }; + let best = (0..5).min_by_key(|&k| cost(&scratch[k])).unwrap_or(0); + out.push(best as u8); + out.extend_from_slice(&scratch[best]); +} + +fn chunk(out: &mut Vec, tag: &[u8; 4], data: &[u8]) { + out.extend_from_slice(&(data.len() as u32).to_be_bytes()); + let start = out.len(); + out.extend_from_slice(tag); + out.extend_from_slice(data); + out.extend_from_slice(&crc32fast::hash(&out[start..]).to_be_bytes()); +} + +/// A PNG of `h` rows of `row_bytes` big-endian sample bytes each (`row(y)` gives row y), with +/// `channels` samples of `depth` bits per pixel, deflated at zlib `level` (0 to 9). +pub fn encode( + h: usize, + w: usize, + channels: usize, + depth: u8, + level: u8, + row: impl Fn(usize, &mut Vec) + Sync, +) -> Vec { + let bpp = channels * depth as usize / 8; + let row_bytes = w * bpp; + let flags = create_comp_flags_from_zip_params(level.min(9) as i32, 0, 0); + let bands: Vec<(Vec, u32, usize)> = (0..h.div_ceil(BAND_ROWS)) + .into_par_iter() + .map(|band| { + let (y0, y1) = (band * BAND_ROWS, ((band + 1) * BAND_ROWS).min(h)); + let mut raw = Vec::with_capacity((y1 - y0) * (row_bytes + 1)); + let mut scratch: [Vec; 5] = std::array::from_fn(|_| Vec::with_capacity(row_bytes)); + let (mut prev, mut cur) = + (Vec::with_capacity(row_bytes), Vec::with_capacity(row_bytes)); + if y0 > 0 { + row(y0 - 1, &mut prev); + } + for y in y0..y1 { + cur.clear(); + row(y, &mut cur); + filter_row( + &cur, + (y > 0).then_some(prev.as_slice()), + bpp, + &mut raw, + &mut scratch, + ); + std::mem::swap(&mut prev, &mut cur); + } + let adler = adler2::adler32_slice(&raw); + let mut c = CompressorOxide::new(flags); + let mut packed = Vec::with_capacity(raw.len() / 2); + let flush = if y1 == h { + TDEFLFlush::Finish + } else { + TDEFLFlush::Sync + }; + let (status, _) = compress_to_output(&mut c, &raw, flush, |o| { + packed.extend_from_slice(o); + true + }); + debug_assert!(matches!(status, TDEFLStatus::Okay | TDEFLStatus::Done)); + (packed, adler, raw.len()) + }) + .collect(); + + let mut png = Vec::with_capacity(bands.iter().map(|b| b.0.len() + 12).sum::() + 64); + png.extend_from_slice(b"\x89PNG\r\n\x1a\n"); + let color = match channels { + 1 => 0u8, + 2 => 4, + 3 => 2, + _ => 6, + }; + let mut ihdr = Vec::with_capacity(13); + ihdr.extend_from_slice(&(w as u32).to_be_bytes()); + ihdr.extend_from_slice(&(h as u32).to_be_bytes()); + ihdr.extend_from_slice(&[depth, color, 0, 0, 0]); + chunk(&mut png, b"IHDR", &ihdr); + let flevel = match level { + 0 | 1 => 0x01u8, + 2..=5 => 0x5e, + 6 => 0x9c, + _ => 0xda, + }; + chunk(&mut png, b"IDAT", &[0x78, flevel]); + let mut adler = 1u32; + for (packed, a, len) in &bands { + chunk(&mut png, b"IDAT", packed); + adler = adler32_combine(adler, *a, *len); + } + chunk(&mut png, b"IDAT", &adler.to_be_bytes()); + chunk(&mut png, b"IEND", &[]); + png +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn adler_combine_matches_one_pass() { + let data: Vec = (0..100_000u32).map(|i| (i * 7 % 251) as u8).collect(); + let (a, b) = data.split_at(31_337); + let whole = adler2::adler32_slice(&data); + let combined = adler32_combine(adler2::adler32_slice(a), adler2::adler32_slice(b), b.len()); + assert_eq!(whole, combined); + } + + #[test] + fn paeth_picks_the_nearest() { + assert_eq!(paeth(10, 20, 10), 20); + assert_eq!(paeth(20, 10, 10), 20); + assert_eq!(paeth(5, 5, 5), 5); + } +} diff --git a/negpy_core/src/py/jpeg.rs b/negpy_core/src/py/jpeg.rs index 2d4a0beee..bbeff37c2 100644 --- a/negpy_core/src/py/jpeg.rs +++ b/negpy_core/src/py/jpeg.rs @@ -41,6 +41,53 @@ fn encode_jpeg<'py>( Ok(PyBytes::new(py, &out)) } +fn dims(shape: &[usize]) -> PyResult<(usize, usize, usize)> { + match shape[..] { + [h, w] => Ok((h, w, 1)), + [h, w, c] if (1..=4).contains(&c) => Ok((h, w, c)), + _ => Err(PyValueError::new_err( + "encode_png needs an HxW or HxWxC array, C 1 to 4", + )), + } +} + +/// A PNG of an HxW or HxWxC (C 1 to 4) uint8 or uint16 array, zlib `level` 0 to 9, filtered +/// and deflated on every core. No ancillary chunks. +#[pyfunction] +#[pyo3(signature = (img, level=6))] +fn encode_png<'py>( + py: Python<'py>, + img: &Bound<'py, PyAny>, + level: u8, +) -> PyResult> { + let png = if let Ok(a) = img.extract::>() { + let (h, w, cn) = dims(a.shape())?; + let src = a + .as_slice() + .map_err(|_| PyValueError::new_err("needs a contiguous array"))?; + super::detach(py, src.len(), || { + crate::png::encode(h, w, cn, 16, level, |y, out| { + for v in &src[y * w * cn..(y + 1) * w * cn] { + out.extend_from_slice(&v.to_be_bytes()); + } + }) + }) + } else { + let a = img.extract::>()?; + let (h, w, cn) = dims(a.shape())?; + let src = a + .as_slice() + .map_err(|_| PyValueError::new_err("needs a contiguous array"))?; + super::detach(py, src.len(), || { + crate::png::encode(h, w, cn, 8, level, |y, out| { + out.extend_from_slice(&src[y * w * cn..(y + 1) * w * cn]) + }) + }) + }; + Ok(PyBytes::new(py, &png)) +} + pub fn register(m: &Bound<'_, PyModule>) -> PyResult<()> { + m.add_function(wrap_pyfunction!(encode_png, m)?)?; m.add_function(wrap_pyfunction!(encode_jpeg, m)?) } diff --git a/tests/test_core_png.py b/tests/test_core_png.py new file mode 100644 index 000000000..68be5e19f --- /dev/null +++ b/tests/test_core_png.py @@ -0,0 +1,35 @@ +import io + +import imagecodecs +import negpy_core +import numpy as np +import pytest +from PIL import Image + + +@pytest.mark.parametrize("shape", [(1, 1), (1, 1, 3), (7, 5, 3), (130, 97, 3), (200, 31), (65, 3, 4), (129, 17, 2)]) +@pytest.mark.parametrize("dtype", [np.uint8, np.uint16]) +@pytest.mark.parametrize("level", [0, 1, 6, 9]) +def test_png_round_trips_through_libpng(shape, dtype, level): + rng = np.random.default_rng(sum(shape) + level) + arr = rng.integers(0, np.iinfo(dtype).max + 1, size=shape, dtype=dtype) + arr[: shape[0] // 2] //= 7 # smooth rows the filters and the deflate can win on + data = bytes(negpy_core.encode_png(arr, level)) + assert np.array_equal(imagecodecs.png_decode(data), arr) + + +def test_png_reads_back_in_pil(): + arr = np.random.default_rng(1).integers(0, 256, size=(300, 211, 3), dtype=np.uint8) + img = Image.open(io.BytesIO(bytes(negpy_core.encode_png(arr, 6)))) + img.load() + assert np.array_equal(np.asarray(img), arr) + + +def test_png_level_zero_is_larger(): + arr = np.tile(np.arange(256, dtype=np.uint8), (300, 2)) + assert len(negpy_core.encode_png(arr, 0)) > len(negpy_core.encode_png(arr, 9)) + + +def test_png_rejects_other_shapes(): + with pytest.raises(ValueError): + negpy_core.encode_png(np.zeros((2, 2, 5), np.uint8), 6) From de15d507838b59ce92771c958221edef63ed9b04 Mon Sep 17 00:00:00 2001 From: Marcin Zawalski Date: Sun, 25 Oct 2026 10:00:51 +0100 Subject: [PATCH 11/12] perf(loaders): keep linear 16-bit scans in uint16 through decode A linear 16-bit TIFF (and Pakon raw) was converted to float32 in the loader, back to uint16 by NonStandardFileWrapper.postprocess in NumPy, and to float32 again by the caller. uint16 -> float32 -> uint16 restores every one of the 65536 codes, so the wrapper now holds the uint16 source, makes float32 only on demand, and a linear 16-bit postprocess() returns the source as it is. Bits are unchanged; tagged sRGB/Adobe and 8-bit sources keep the float path. 23.5 MP TIFF, load_linear_preview without IR: ~330 -> ~295 ms; every export decode of such a scan skips the same two passes. --- negpy/infrastructure/loaders/helpers.py | 25 ++++++++++++++++---- negpy/infrastructure/loaders/pakon_loader.py | 3 +-- negpy/infrastructure/loaders/tiff_loader.py | 20 +++++++++------- 3 files changed, 34 insertions(+), 14 deletions(-) diff --git a/negpy/infrastructure/loaders/helpers.py b/negpy/infrastructure/loaders/helpers.py index 404b83d34..3a186a012 100644 --- a/negpy/infrastructure/loaders/helpers.py +++ b/negpy/infrastructure/loaders/helpers.py @@ -12,7 +12,7 @@ from negpy.domain.models import ColorSpace from negpy.features.process.models import DemosaicMode from negpy.infrastructure.loaders.constants import SUPPORTED_RAW_EXTENSIONS -from negpy.kernel.image.logic import apply_exif_orientation, ensure_rgb +from negpy.kernel.image.logic import apply_exif_orientation, ensure_rgb, uint16_to_float32 from negpy.kernel.system.config import APP_CONFIG from negpy.kernel.system.logging import get_logger @@ -649,11 +649,16 @@ class NonStandardFileWrapper: def __init__( self, - data: np.ndarray, + data: Optional[np.ndarray] = None, full_output_hw: Optional[Tuple[int, int]] = None, wb_gains: Optional[Tuple[float, float, float]] = None, + *, + source_u16: Optional[np.ndarray] = None, ) -> None: - self.data = data + self._data = data + # A linear 16-bit source: float32 is made only on demand, and a linear 16-bit + # postprocess() hands it back as it is. uint16 -> float32 -> uint16 restores every code. + self._source_u16 = source_u16 # If set, `sizes` reports this (h, w) for full image; else derived from `data` shape. self._full_output_hw: Optional[Tuple[int, int]] = full_output_hw # As-shot (R, G, B) white balance gains, applied by postprocess() when the caller asks @@ -661,12 +666,21 @@ def __init__( # DNGs, which are always neutral. postprocess() then leaves the data untouched. self.wb_gains: Optional[Tuple[float, float, float]] = wb_gains + @property + def data(self) -> np.ndarray: + if self._data is None: + assert self._source_u16 is not None + self._data = uint16_to_float32(self._source_u16) + return self._data + @property def sizes(self) -> Any: if self._full_output_hw is not None: h, w = self._full_output_hw else: - h, w = self.data.shape[0], self.data.shape[1] + source = self._data if self._data is not None else self._source_u16 + assert source is not None + h, w = source.shape[0], source.shape[1] return SimpleNamespace(raw_height=int(h), raw_width=int(w)) def __enter__(self) -> "NonStandardFileWrapper": @@ -679,6 +693,9 @@ def postprocess(self, **kwargs: Any) -> np.ndarray: bps = kwargs.get("output_bps", 8) half_size = kwargs.get("half_size", False) gamma = kwargs.get("gamma") + camera_wb = kwargs.get("use_camera_wb") and self.wb_gains is not None + if self._source_u16 is not None and bps == 16 and gamma is not None and tuple(gamma) == (1, 1) and not camera_wb: + return np.ascontiguousarray(self._source_u16[::2, ::2] if half_size else self._source_u16) data = self.data if half_size: data = data[::2, ::2] diff --git a/negpy/infrastructure/loaders/pakon_loader.py b/negpy/infrastructure/loaders/pakon_loader.py index 112e47bdb..a55cdf516 100644 --- a/negpy/infrastructure/loaders/pakon_loader.py +++ b/negpy/infrastructure/loaders/pakon_loader.py @@ -6,7 +6,6 @@ from PIL import Image from negpy.domain.interfaces import IImageLoader from negpy.infrastructure.loaders.helpers import NonStandardFileWrapper, linear_uint16_to_display_uint8 -from negpy.kernel.image.logic import uint16_to_float32 class PakonLoader(IImageLoader): @@ -97,7 +96,7 @@ def load(self, file_path: str) -> Tuple[ContextManager[Any], dict]: data = data.reshape((h, w, 3))[..., ::-1] metadata = {"orientation": 0, "ir": None} - return NonStandardFileWrapper(uint16_to_float32(np.ascontiguousarray(data))), metadata + return NonStandardFileWrapper(source_u16=np.ascontiguousarray(data)), metadata except Exception as e: # Fallback to Rawpy or re-raise to be caught by worker raise RuntimeError(f"Pakon Load Failure: {e}") from e diff --git a/negpy/infrastructure/loaders/tiff_loader.py b/negpy/infrastructure/loaders/tiff_loader.py index c00bff386..2ab23c9b9 100644 --- a/negpy/infrastructure/loaders/tiff_loader.py +++ b/negpy/infrastructure/loaders/tiff_loader.py @@ -203,13 +203,6 @@ def load( if ir is None and load_ir: ir = _read_ir_from_extra_page(file_path, img.shape[0], img.shape[1]) - if img.dtype == np.uint8: - f32 = uint8_to_float32(np.ascontiguousarray(img)) - elif img.dtype == np.uint16: - f32 = uint16_to_float32(np.ascontiguousarray(img)) - else: - f32 = np.clip(img.astype(np.float32), 0, 1) - icc_bytes: bytes | None = None try: with tifffile.TiffFile(file_path) as tif: @@ -229,12 +222,23 @@ def load( # linear, which no ColorSpace names, so it stays None; a positive source has # resolved that ambiguity and takes the 8-bit assumption. color_space = ColorSpace.SRGB.value + decodes = color_space in (ColorSpace.SRGB.value, ColorSpace.ADOBE_RGB.value) + if img.dtype == np.uint16 and not decodes: + wrapper = NonStandardFileWrapper(source_u16=np.ascontiguousarray(img)) + else: + if img.dtype == np.uint8: + f32 = uint8_to_float32(np.ascontiguousarray(img)) + elif img.dtype == np.uint16: + f32 = uint16_to_float32(np.ascontiguousarray(img)) + else: + f32 = np.clip(img.astype(np.float32), 0, 1) if color_space == ColorSpace.SRGB.value: f32 = srgb_to_linear(f32) elif color_space == ColorSpace.ADOBE_RGB.value: # Adobe RGB's TRC is the working space's own gamma, so its decode is the # inverse of the pipeline OETF encode. f32 = working_oetf_decode(f32) + wrapper = NonStandardFileWrapper(f32) metadata = { "orientation": read_orientation(file_path), "color_space": color_space, @@ -242,7 +246,7 @@ def load( "ir": ir, "ir_valid_mask": ir_valid_mask, } - return NonStandardFileWrapper(f32), metadata + return wrapper, metadata def load_bounded_preview( self, From 2c8d97a8a6d8b557d02c2cc0cb203e262205a648 Mon Sep 17 00:00:00 2001 From: Marcin Zawalski Date: Sun, 25 Oct 2026 13:28:32 +0100 Subject: [PATCH 12/12] fix(retouch): dust density log in Rust, same bits on every CPU test_compute_dust_stats_pin failed on some CI runners. The density plane used NumPy's float32 log10, whose SIMD code differs between AVX-512 and AVX2 CPUs, so the pinned output depended on the runner. The density now comes from negpy_core.log_density_open (correctly rounded log10f). The new pin equals NumPy's AVX2 result and holds with NumPy's AVX-512 paths on or off. --- negpy/features/retouch/logic.py | 2 +- tests/test_core_ir.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/negpy/features/retouch/logic.py b/negpy/features/retouch/logic.py index cd185df9b..2f7dd24d3 100644 --- a/negpy/features/retouch/logic.py +++ b/negpy/features/retouch/logic.py @@ -225,7 +225,7 @@ def _density(img: ImageBuffer) -> np.ndarray: """Log density of the luminance. Local contrast measured here is exposure-invariant: a gain on the source is an offset in density and cancels in every difference.""" - return -np.log10(np.clip(get_luminance(img), 1e-6, None)).astype(np.float32) + return -negpy_core.log_density_open(get_luminance(img)) def _mask_to_score(mask: np.ndarray, pad_px: float) -> np.ndarray: diff --git a/tests/test_core_ir.py b/tests/test_core_ir.py index fd83bc592..60f5f047d 100644 --- a/tests/test_core_ir.py +++ b/tests/test_core_ir.py @@ -102,7 +102,7 @@ def test_compute_dust_stats_pin(): rng = np.random.default_rng(8) img = (0.3 + 0.2 * rng.random((140, 180, 3))).astype(np.float32) img[40:43, 50:53] *= 0.2 # a speck - assert digest(L.compute_dust_stats(img, 2)) == "1e3e0ae2a35da752b5c5e36fa2286c8e59d1f00d9e6aa26e504468e89c85ecdc" + assert digest(L.compute_dust_stats(img, 2)) == "d0ebd32203ccb8dc9e9197e444323b30be9bcf2ebeee3565b2f4d05762b3e320" def test_split_hairs_matches_components():