From e21537c409c1237ebe189aaed83c8d6fb5444b41 Mon Sep 17 00:00:00 2001 From: Constantin Pape Date: Sun, 17 May 2026 15:21:55 +0200 Subject: [PATCH 1/3] First version of filters --- CMakeLists.txt | 1 + MIGRATION_GUIDE.md | 82 +++ include/bioimage_cpp/filters/convolve.hxx | 451 +++++++++++++++ include/bioimage_cpp/filters/eigenvalues.hxx | 134 +++++ include/bioimage_cpp/filters/gaussian.hxx | 494 ++++++++++++++++ include/bioimage_cpp/filters/kernel.hxx | 137 +++++ src/bindings/filters.cxx | 565 +++++++++++++++++++ src/bindings/filters.hxx | 9 + src/bindings/module.cxx | 2 + src/bioimage_cpp/__init__.py | 2 + src/bioimage_cpp/filters/__init__.py | 20 + src/bioimage_cpp/filters/_filters.py | 236 ++++++++ tests/test_filters.py | 271 +++++++++ 13 files changed, 2404 insertions(+) create mode 100644 include/bioimage_cpp/filters/convolve.hxx create mode 100644 include/bioimage_cpp/filters/eigenvalues.hxx create mode 100644 include/bioimage_cpp/filters/gaussian.hxx create mode 100644 include/bioimage_cpp/filters/kernel.hxx create mode 100644 src/bindings/filters.cxx create mode 100644 src/bindings/filters.hxx create mode 100644 src/bioimage_cpp/filters/__init__.py create mode 100644 src/bioimage_cpp/filters/_filters.py create mode 100644 tests/test_filters.py diff --git a/CMakeLists.txt b/CMakeLists.txt index e56cde5..e75d985 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -13,6 +13,7 @@ nanobind_add_module(_core NB_STATIC src/bindings/blocking.cxx src/bindings/module.cxx + src/bindings/filters.cxx src/bindings/graph.cxx src/bindings/ground_truth.cxx src/bindings/segmentation.cxx diff --git a/MIGRATION_GUIDE.md b/MIGRATION_GUIDE.md index 6ac442b..9777b1d 100644 --- a/MIGRATION_GUIDE.md +++ b/MIGRATION_GUIDE.md @@ -946,6 +946,88 @@ Notes: - Every value in the input must be present in the mapping. - Non-contiguous inputs are copied before entering C++. +## Image Filters + +`bioimage-cpp` ships a small Gaussian-derivative filter set under +`bic.filters`. The scope is the "ilastik filter set" exposed by +`fastfilters`, which is also the most-used subset of `vigra.filters`. + +Vigra / fastfilters: + +```python +import vigra.filters as vf +out = vf.gaussianSmoothing(img, sigma=1.5) +ev = vf.hessianOfGaussianEigenvalues(img, scale=1.5) + +import fastfilters as ff +out = ff.gaussianSmoothing(img, sigma=1.5) +ev = ff.hessianOfGaussianEigenvalues(img, scale=1.5) +``` + +bioimage-cpp: + +```python +import bioimage_cpp as bic + +out = bic.filters.gaussian_smoothing(img, sigma=1.5) +ev = bic.filters.hessian_of_gaussian_eigenvalues(img, sigma=1.5) +``` + +Name mapping: + +| vigra / fastfilters name | bioimage-cpp name | +| --- | --- | +| `gaussianSmoothing` | `gaussian_smoothing` | +| `gaussianDerivative` | `gaussian_derivative` | +| `gaussianGradientMagnitude` | `gaussian_gradient_magnitude` | +| `laplacianOfGaussian` | `laplacian_of_gaussian` | +| `hessianOfGaussianEigenvalues` | `hessian_of_gaussian_eigenvalues` | +| `structureTensorEigenvalues` | `structure_tensor_eigenvalues` | + +Common parameters: + +- `sigma` is a positive scalar or a per-axis sequence of length + `image.ndim`. Anisotropic sigma is supported on every filter. +- `gaussian_derivative` takes an `order` argument that is a scalar or a + per-axis sequence of ints in `{0, 1, 2}`. +- `structure_tensor_eigenvalues` takes positional `inner_sigma` and + `outer_sigma` (vigra calls them `innerScale` / `outerScale`). +- `window_size` controls the kernel radius: + `radius = ceil(window_size * sigma)`. `0.0` (the default) selects the + vigra-style default `3 + 0.5 * order`. Matches the same-named parameter + in vigra/fastfilters. + +Important differences from vigra and fastfilters: + +- Only 2D and 3D scalar (single-channel) inputs are supported in v1. + Channels and leading batch axes should be looped externally — matches + fastfilters' convention. Vigra's `taggedView`/`AxisInfo` machinery is + not reproduced. +- C++ kernels operate on `float32`. `float64` inputs are accepted and the + output is cast back to `float64`. `uint8` and `uint16` are accepted with + a `float32` output (the typical ML-feature use case). +- Boundary handling is `mirror` (matches scipy `mode="mirror"` — + reflection without edge-pixel repeat). Other boundary modes are not + exposed yet; the C++ layer carries an enum for future tiled processing. +- Eigenvalue outputs have a trailing axis of size `image.ndim`, sorted + largest → smallest. This matches `fastfilters`. To get vigra's + ascending order, reverse with `result[..., ::-1]`. +- No IIR / recursive Gaussian, no `convolve` / `recursiveFilter2D`, no + morphology, no distance transforms, no nonlinear diffusion, and no + non-local means in v1. Use `scipy.ndimage`, `skimage`, or the original + vigra/fastfilters bindings if you need those. + +Implementation notes: + +- All six filters are written as portable C++20 scalar code that the + compiler auto-vectorizes. No SIMD intrinsics, no per-file ISA flags, no + runtime CPU dispatch, no vendored SIMD library. This keeps the build + light enough to ship as portable PyPI wheels across Linux/macOS/Windows + and x86_64/arm64. +- Single-threaded for now. Threading can be added later via + `detail/threading.hxx::parallel_for_chunks` without changing the + public API. + ## I/O and Build Dependencies `bioimage-cpp` intentionally does not replace nifty or affogato I/O helpers. diff --git a/include/bioimage_cpp/filters/convolve.hxx b/include/bioimage_cpp/filters/convolve.hxx new file mode 100644 index 0000000..c25a4ee --- /dev/null +++ b/include/bioimage_cpp/filters/convolve.hxx @@ -0,0 +1,451 @@ +#pragma once + +#include "bioimage_cpp/filters/kernel.hxx" + +#include +#include + +namespace bioimage_cpp::filters { + +namespace detail { + +// Mirror reflection without edge-pixel repeat (scipy mode="mirror"). +// Period is 2*(n-1). Handles arbitrary integer x. +inline std::ptrdiff_t mirror_index(std::ptrdiff_t x, std::ptrdiff_t n) { + if (n <= 1) return 0; + const std::ptrdiff_t period = 2 * (n - 1); + std::ptrdiff_t r = x % period; + if (r < 0) r += period; + if (r >= n) r = period - r; + return r; +} + +// Convolve along the contiguous (innermost) axis. Specialised for compile-time +// radius R and symmetry. +template +void convolve_x_radius( + const float *__restrict in, + float *__restrict out, + std::ptrdiff_t n_rows, + std::ptrdiff_t n_cols, + const float *__restrict h +) { + if (n_cols <= 0 || n_rows <= 0) return; + + const std::ptrdiff_t prologue_end = std::min(R, n_cols); + const std::ptrdiff_t epilogue_start = std::max(prologue_end, n_cols - R); + + for (std::ptrdiff_t row = 0; row < n_rows; ++row) { + const float *__restrict in_row = in + row * n_cols; + float *__restrict out_row = out + row * n_cols; + + // Border prologue + for (std::ptrdiff_t x = 0; x < prologue_end; ++x) { + float acc; + if constexpr (Symmetric) { + acc = h[0] * in_row[x]; + for (int k = 1; k <= R; ++k) { + const float left = in_row[mirror_index(x - k, n_cols)]; + const float right = in_row[mirror_index(x + k, n_cols)]; + acc += h[k] * (left + right); + } + } else { + acc = 0.0f; + for (int k = 1; k <= R; ++k) { + const float left = in_row[mirror_index(x - k, n_cols)]; + const float right = in_row[mirror_index(x + k, n_cols)]; + acc += h[k] * (right - left); + } + } + out_row[x] = acc; + } + + // Main loop (branchless): the compiler vectorizes the outer x loop. + if constexpr (Symmetric) { + for (std::ptrdiff_t x = prologue_end; x < epilogue_start; ++x) { + float acc = h[0] * in_row[x]; + for (int k = 1; k <= R; ++k) { + acc += h[k] * (in_row[x + k] + in_row[x - k]); + } + out_row[x] = acc; + } + } else { + for (std::ptrdiff_t x = prologue_end; x < epilogue_start; ++x) { + float acc = 0.0f; + for (int k = 1; k <= R; ++k) { + acc += h[k] * (in_row[x + k] - in_row[x - k]); + } + out_row[x] = acc; + } + } + + // Border epilogue + for (std::ptrdiff_t x = epilogue_start; x < n_cols; ++x) { + float acc; + if constexpr (Symmetric) { + acc = h[0] * in_row[x]; + for (int k = 1; k <= R; ++k) { + const float left = in_row[mirror_index(x - k, n_cols)]; + const float right = in_row[mirror_index(x + k, n_cols)]; + acc += h[k] * (left + right); + } + } else { + acc = 0.0f; + for (int k = 1; k <= R; ++k) { + const float left = in_row[mirror_index(x - k, n_cols)]; + const float right = in_row[mirror_index(x + k, n_cols)]; + acc += h[k] * (right - left); + } + } + out_row[x] = acc; + } + } +} + +// Same as above but with runtime radius (handles radii > the template cap). +template +void convolve_x_runtime( + const float *__restrict in, + float *__restrict out, + std::ptrdiff_t n_rows, + std::ptrdiff_t n_cols, + int radius, + const float *__restrict h +) { + if (n_cols <= 0 || n_rows <= 0) return; + + const std::ptrdiff_t R = radius; + const std::ptrdiff_t prologue_end = std::min(R, n_cols); + const std::ptrdiff_t epilogue_start = std::max(prologue_end, n_cols - R); + + for (std::ptrdiff_t row = 0; row < n_rows; ++row) { + const float *__restrict in_row = in + row * n_cols; + float *__restrict out_row = out + row * n_cols; + + for (std::ptrdiff_t x = 0; x < prologue_end; ++x) { + float acc; + if constexpr (Symmetric) { + acc = h[0] * in_row[x]; + for (int k = 1; k <= R; ++k) { + acc += h[k] * (in_row[mirror_index(x - k, n_cols)] + + in_row[mirror_index(x + k, n_cols)]); + } + } else { + acc = 0.0f; + for (int k = 1; k <= R; ++k) { + acc += h[k] * (in_row[mirror_index(x + k, n_cols)] - + in_row[mirror_index(x - k, n_cols)]); + } + } + out_row[x] = acc; + } + + if constexpr (Symmetric) { + for (std::ptrdiff_t x = prologue_end; x < epilogue_start; ++x) { + float acc = h[0] * in_row[x]; + for (int k = 1; k <= R; ++k) { + acc += h[k] * (in_row[x + k] + in_row[x - k]); + } + out_row[x] = acc; + } + } else { + for (std::ptrdiff_t x = prologue_end; x < epilogue_start; ++x) { + float acc = 0.0f; + for (int k = 1; k <= R; ++k) { + acc += h[k] * (in_row[x + k] - in_row[x - k]); + } + out_row[x] = acc; + } + } + + for (std::ptrdiff_t x = epilogue_start; x < n_cols; ++x) { + float acc; + if constexpr (Symmetric) { + acc = h[0] * in_row[x]; + for (int k = 1; k <= R; ++k) { + acc += h[k] * (in_row[mirror_index(x - k, n_cols)] + + in_row[mirror_index(x + k, n_cols)]); + } + } else { + acc = 0.0f; + for (int k = 1; k <= R; ++k) { + acc += h[k] * (in_row[mirror_index(x + k, n_cols)] - + in_row[mirror_index(x - k, n_cols)]); + } + } + out_row[x] = acc; + } + } +} + +// X-strip block size for the strided pass. 64 floats = 256 bytes, fits in L1 +// and large enough to amortise the kernel-coefficient broadcast. +inline constexpr std::ptrdiff_t kStripBlock = 64; + +// Convolve along a strided (non-innermost) axis. Logical shape +// (n_outer, n_axis, n_inner) in C-order; we accumulate into strips of +// kStripBlock contiguous columns of the innermost axis, which gives the +// compiler something easy to vectorise. +template +void convolve_strided_radius( + const float *__restrict in, + float *__restrict out, + std::ptrdiff_t n_outer, + std::ptrdiff_t n_axis, + std::ptrdiff_t n_inner, + const float *__restrict h +) { + if (n_axis <= 0 || n_inner <= 0 || n_outer <= 0) return; + + const std::ptrdiff_t outer_stride = n_axis * n_inner; + const std::ptrdiff_t prologue_end = std::min(R, n_axis); + const std::ptrdiff_t epilogue_start = std::max(prologue_end, n_axis - R); + + for (std::ptrdiff_t o = 0; o < n_outer; ++o) { + const float *__restrict in_o = in + o * outer_stride; + float *__restrict out_o = out + o * outer_stride; + + // Main rows (no border on the axis dimension). + for (std::ptrdiff_t y = prologue_end; y < epilogue_start; ++y) { + float *__restrict out_row = out_o + y * n_inner; + const float *__restrict center = in_o + y * n_inner; + + for (std::ptrdiff_t xb = 0; xb < n_inner; xb += kStripBlock) { + const std::ptrdiff_t strip = std::min(kStripBlock, n_inner - xb); + float acc[kStripBlock]; + + if constexpr (Symmetric) { + const float h0 = h[0]; + for (std::ptrdiff_t i = 0; i < strip; ++i) { + acc[i] = h0 * center[xb + i]; + } + } else { + for (std::ptrdiff_t i = 0; i < strip; ++i) { + acc[i] = 0.0f; + } + } + + for (int k = 1; k <= R; ++k) { + const float hk = h[k]; + const float *__restrict up = in_o + (y - k) * n_inner + xb; + const float *__restrict dn = in_o + (y + k) * n_inner + xb; + if constexpr (Symmetric) { + for (std::ptrdiff_t i = 0; i < strip; ++i) { + acc[i] += hk * (up[i] + dn[i]); + } + } else { + for (std::ptrdiff_t i = 0; i < strip; ++i) { + acc[i] += hk * (dn[i] - up[i]); + } + } + } + + for (std::ptrdiff_t i = 0; i < strip; ++i) { + out_row[xb + i] = acc[i]; + } + } + } + + // Border rows (axis prologue + epilogue), with mirror on y±k. + auto run_border_row = [&](std::ptrdiff_t y) { + float *__restrict out_row = out_o + y * n_inner; + const float *__restrict center = in_o + y * n_inner; + + for (std::ptrdiff_t xb = 0; xb < n_inner; xb += kStripBlock) { + const std::ptrdiff_t strip = std::min(kStripBlock, n_inner - xb); + float acc[kStripBlock]; + + if constexpr (Symmetric) { + const float h0 = h[0]; + for (std::ptrdiff_t i = 0; i < strip; ++i) { + acc[i] = h0 * center[xb + i]; + } + } else { + for (std::ptrdiff_t i = 0; i < strip; ++i) { + acc[i] = 0.0f; + } + } + + for (int k = 1; k <= R; ++k) { + const float hk = h[k]; + const std::ptrdiff_t y_up = mirror_index(y - k, n_axis); + const std::ptrdiff_t y_dn = mirror_index(y + k, n_axis); + const float *__restrict up = in_o + y_up * n_inner + xb; + const float *__restrict dn = in_o + y_dn * n_inner + xb; + if constexpr (Symmetric) { + for (std::ptrdiff_t i = 0; i < strip; ++i) { + acc[i] += hk * (up[i] + dn[i]); + } + } else { + for (std::ptrdiff_t i = 0; i < strip; ++i) { + acc[i] += hk * (dn[i] - up[i]); + } + } + } + + for (std::ptrdiff_t i = 0; i < strip; ++i) { + out_row[xb + i] = acc[i]; + } + } + }; + + for (std::ptrdiff_t y = 0; y < prologue_end; ++y) run_border_row(y); + for (std::ptrdiff_t y = epilogue_start; y < n_axis; ++y) run_border_row(y); + } +} + +// Runtime-radius fallback for the strided pass. Falls back to per-pixel +// indexing rather than the strip pattern; intended only for the rare case +// where the kernel radius exceeds the compile-time cap. +template +void convolve_strided_runtime( + const float *__restrict in, + float *__restrict out, + std::ptrdiff_t n_outer, + std::ptrdiff_t n_axis, + std::ptrdiff_t n_inner, + int radius, + const float *__restrict h +) { + if (n_axis <= 0 || n_inner <= 0 || n_outer <= 0) return; + + const std::ptrdiff_t R = radius; + const std::ptrdiff_t outer_stride = n_axis * n_inner; + + for (std::ptrdiff_t o = 0; o < n_outer; ++o) { + const float *__restrict in_o = in + o * outer_stride; + float *__restrict out_o = out + o * outer_stride; + + for (std::ptrdiff_t y = 0; y < n_axis; ++y) { + float *__restrict out_row = out_o + y * n_inner; + const float *__restrict center = in_o + y * n_inner; + const bool needs_mirror = (y < R) || (y >= n_axis - R); + + for (std::ptrdiff_t xb = 0; xb < n_inner; xb += kStripBlock) { + const std::ptrdiff_t strip = std::min(kStripBlock, n_inner - xb); + float acc[kStripBlock]; + + if constexpr (Symmetric) { + const float h0 = h[0]; + for (std::ptrdiff_t i = 0; i < strip; ++i) acc[i] = h0 * center[xb + i]; + } else { + for (std::ptrdiff_t i = 0; i < strip; ++i) acc[i] = 0.0f; + } + + for (int k = 1; k <= R; ++k) { + const float hk = h[k]; + const std::ptrdiff_t y_up = + needs_mirror ? mirror_index(y - k, n_axis) : (y - k); + const std::ptrdiff_t y_dn = + needs_mirror ? mirror_index(y + k, n_axis) : (y + k); + const float *up = in_o + y_up * n_inner + xb; + const float *dn = in_o + y_dn * n_inner + xb; + if constexpr (Symmetric) { + for (std::ptrdiff_t i = 0; i < strip; ++i) { + acc[i] += hk * (up[i] + dn[i]); + } + } else { + for (std::ptrdiff_t i = 0; i < strip; ++i) { + acc[i] += hk * (dn[i] - up[i]); + } + } + } + + for (std::ptrdiff_t i = 0; i < strip; ++i) out_row[xb + i] = acc[i]; + } + } + } +} + +} // namespace detail + +// Maximum compile-time-specialised kernel radius. Kernels with radius > this +// dispatch to the runtime-radius fallback. 12 covers sigma up to ~3.4 with the +// default window (3*sigma); larger sigma is supported but slower. +inline constexpr int kMaxSpecialisedRadius = 12; + +// Convolve along the innermost (contiguous) axis. in and out must not alias. +inline void convolve_axis_x( + const float *in, + float *out, + std::ptrdiff_t n_rows, + std::ptrdiff_t n_cols, + const Kernel1D &kernel +) { + const int r = kernel.radius; + const float *h = kernel.half_coefs.data(); + const bool sym = kernel.is_symmetric; + +#define BIOIMAGE_FILTERS_DISPATCH_RADIUS_X(R) \ + case R: \ + if (sym) detail::convolve_x_radius(in, out, n_rows, n_cols, h); \ + else detail::convolve_x_radius(in, out, n_rows, n_cols, h); \ + return; + + switch (r) { + BIOIMAGE_FILTERS_DISPATCH_RADIUS_X(1) + BIOIMAGE_FILTERS_DISPATCH_RADIUS_X(2) + BIOIMAGE_FILTERS_DISPATCH_RADIUS_X(3) + BIOIMAGE_FILTERS_DISPATCH_RADIUS_X(4) + BIOIMAGE_FILTERS_DISPATCH_RADIUS_X(5) + BIOIMAGE_FILTERS_DISPATCH_RADIUS_X(6) + BIOIMAGE_FILTERS_DISPATCH_RADIUS_X(7) + BIOIMAGE_FILTERS_DISPATCH_RADIUS_X(8) + BIOIMAGE_FILTERS_DISPATCH_RADIUS_X(9) + BIOIMAGE_FILTERS_DISPATCH_RADIUS_X(10) + BIOIMAGE_FILTERS_DISPATCH_RADIUS_X(11) + BIOIMAGE_FILTERS_DISPATCH_RADIUS_X(12) + default: + if (sym) detail::convolve_x_runtime(in, out, n_rows, n_cols, r, h); + else detail::convolve_x_runtime(in, out, n_rows, n_cols, r, h); + return; + } +#undef BIOIMAGE_FILTERS_DISPATCH_RADIUS_X +} + +// Convolve along a strided axis. Logical layout is (n_outer, n_axis, n_inner) +// in C-order; in and out must not alias. +inline void convolve_axis_strided( + const float *in, + float *out, + std::ptrdiff_t n_outer, + std::ptrdiff_t n_axis, + std::ptrdiff_t n_inner, + const Kernel1D &kernel +) { + const int r = kernel.radius; + const float *h = kernel.half_coefs.data(); + const bool sym = kernel.is_symmetric; + +#define BIOIMAGE_FILTERS_DISPATCH_RADIUS_S(R) \ + case R: \ + if (sym) \ + detail::convolve_strided_radius(in, out, n_outer, n_axis, n_inner, h); \ + else \ + detail::convolve_strided_radius(in, out, n_outer, n_axis, n_inner, h); \ + return; + + switch (r) { + BIOIMAGE_FILTERS_DISPATCH_RADIUS_S(1) + BIOIMAGE_FILTERS_DISPATCH_RADIUS_S(2) + BIOIMAGE_FILTERS_DISPATCH_RADIUS_S(3) + BIOIMAGE_FILTERS_DISPATCH_RADIUS_S(4) + BIOIMAGE_FILTERS_DISPATCH_RADIUS_S(5) + BIOIMAGE_FILTERS_DISPATCH_RADIUS_S(6) + BIOIMAGE_FILTERS_DISPATCH_RADIUS_S(7) + BIOIMAGE_FILTERS_DISPATCH_RADIUS_S(8) + BIOIMAGE_FILTERS_DISPATCH_RADIUS_S(9) + BIOIMAGE_FILTERS_DISPATCH_RADIUS_S(10) + BIOIMAGE_FILTERS_DISPATCH_RADIUS_S(11) + BIOIMAGE_FILTERS_DISPATCH_RADIUS_S(12) + default: + if (sym) + detail::convolve_strided_runtime(in, out, n_outer, n_axis, n_inner, r, h); + else + detail::convolve_strided_runtime(in, out, n_outer, n_axis, n_inner, r, h); + return; + } +#undef BIOIMAGE_FILTERS_DISPATCH_RADIUS_S +} + +} // namespace bioimage_cpp::filters diff --git a/include/bioimage_cpp/filters/eigenvalues.hxx b/include/bioimage_cpp/filters/eigenvalues.hxx new file mode 100644 index 0000000..0bece2f --- /dev/null +++ b/include/bioimage_cpp/filters/eigenvalues.hxx @@ -0,0 +1,134 @@ +#pragma once + +#include +#include +#include + +namespace bioimage_cpp::filters { + +// Eigenvalues of a 2x2 symmetric matrix +// [xx xy] +// [xy yy] +// computed via the closed form (tr/2) +- sqrt(((xx-yy)/2)^2 + xy^2), then +// sorted descending. SoA inputs/outputs of length n. +inline void ev2_symmetric_descending( + const float *__restrict xx, + const float *__restrict xy, + const float *__restrict yy, + float *__restrict e_large, + float *__restrict e_small, + std::ptrdiff_t n +) { + for (std::ptrdiff_t i = 0; i < n; ++i) { + const float half_tr = 0.5f * (xx[i] + yy[i]); + const float half_diff = 0.5f * (xx[i] - yy[i]); + const float disc = std::sqrt(half_diff * half_diff + xy[i] * xy[i]); + e_large[i] = half_tr + disc; + e_small[i] = half_tr - disc; + } +} + +namespace detail { + +// Eigenvalues of a 3x3 symmetric matrix at a single pixel, sorted descending. +// Uses David Eberly's trigonometric closed-form solution (see VIGRA's +// mathutil.hxx). Steps: normalise by the maximum absolute element to guard +// against overflow; compute the deviatoric matrix B = A - (tr/3) I; the three +// eigenvalues are a + 2*p*cos(phi + k*2pi/3) where p = sqrt(tr(B^2)/6), +// phi = acos(clamp(det(B)/(2 p^3), -1, 1))/3. +inline void ev3_one_descending( + float a00, float a01, float a02, + float a11, float a12, float a22, + float &e0, float &e1, float &e2 +) { + float m = std::abs(a00); + m = std::max(m, std::abs(a01)); + m = std::max(m, std::abs(a02)); + m = std::max(m, std::abs(a11)); + m = std::max(m, std::abs(a12)); + m = std::max(m, std::abs(a22)); + + if (m == 0.0f) { + e0 = 0.0f; + e1 = 0.0f; + e2 = 0.0f; + return; + } + + const float inv_m = 1.0f / m; + a00 *= inv_m; + a01 *= inv_m; + a02 *= inv_m; + a11 *= inv_m; + a12 *= inv_m; + a22 *= inv_m; + + const float a = (a00 + a11 + a22) * (1.0f / 3.0f); + const float b00 = a00 - a; + const float b11 = a11 - a; + const float b22 = a22 - a; + + // trace(B B^T) = b00^2 + b11^2 + b22^2 + 2*(a01^2 + a02^2 + a12^2) + const float trace_b2 = + b00 * b00 + b11 * b11 + b22 * b22 + 2.0f * (a01 * a01 + a02 * a02 + a12 * a12); + const float p2 = trace_b2 * (1.0f / 6.0f); + + if (!(p2 > 0.0f)) { + const float val = a * m; + e0 = val; + e1 = val; + e2 = val; + return; + } + const float p = std::sqrt(p2); + + const float det_b = b00 * (b11 * b22 - a12 * a12) + - a01 * (a01 * b22 - a12 * a02) + + a02 * (a01 * a12 - b11 * a02); + + float r = det_b / (2.0f * p2 * p); + // Clamp to [-1, 1] to guard against floating-point excursions that would + // otherwise produce NaN from acos. + r = std::max(-1.0f, std::min(1.0f, r)); + + constexpr float kTwoPiOver3 = 2.0943951023931953f; + const float phi = std::acos(r) * (1.0f / 3.0f); + const float two_p = 2.0f * p; + + const float root_large = a + two_p * std::cos(phi); + const float root_small = a + two_p * std::cos(phi + kTwoPiOver3); + const float root_mid = 3.0f * a - root_large - root_small; + + e0 = m * root_large; + e1 = m * root_mid; + e2 = m * root_small; +} + +} // namespace detail + +// Eigenvalues of a 3x3 symmetric matrix +// [xx xy xz] +// [xy yy yz] +// [xz yz zz] +// sorted descending (e0 >= e1 >= e2). SoA inputs/outputs of length n. +inline void ev3_symmetric_descending( + const float *__restrict xx, + const float *__restrict xy, + const float *__restrict xz, + const float *__restrict yy, + const float *__restrict yz, + const float *__restrict zz, + float *__restrict e0, + float *__restrict e1, + float *__restrict e2, + std::ptrdiff_t n +) { + for (std::ptrdiff_t i = 0; i < n; ++i) { + detail::ev3_one_descending( + xx[i], xy[i], xz[i], yy[i], yz[i], zz[i], + e0[i], e1[i], e2[i] + ); + } +} + +} // namespace bioimage_cpp::filters diff --git a/include/bioimage_cpp/filters/gaussian.hxx b/include/bioimage_cpp/filters/gaussian.hxx new file mode 100644 index 0000000..0893c7a --- /dev/null +++ b/include/bioimage_cpp/filters/gaussian.hxx @@ -0,0 +1,494 @@ +#pragma once + +#include "bioimage_cpp/filters/convolve.hxx" +#include "bioimage_cpp/filters/eigenvalues.hxx" +#include "bioimage_cpp/filters/kernel.hxx" + +#include +#include +#include +#include +#include + +namespace bioimage_cpp::filters { + +// --------------------------------------------------------------------------- +// Separable Gaussian (derivative) along each axis. +// --------------------------------------------------------------------------- +// +// `in`, `out`, `workspace` are C-contiguous buffers of the same total size. +// `out` and `workspace` must NOT alias `in` (the binding ensures this by +// allocating fresh output and scratch). `out` and `workspace` may not alias +// each other. After the call, `out` holds the filter response and `workspace` +// is left in an unspecified state. + +inline void gaussian_separable_2d( + const float *in, + float *out, + float *workspace, + std::ptrdiff_t ny, + std::ptrdiff_t nx, + const Kernel1D &ky, + const Kernel1D &kx +) { + // Axis 0 (Y), strided: (n_outer=1, n_axis=ny, n_inner=nx). + convolve_axis_strided(in, workspace, 1, ny, nx, ky); + // Axis 1 (X), contiguous: (n_rows=ny, n_cols=nx). + convolve_axis_x(workspace, out, ny, nx, kx); +} + +inline void gaussian_separable_3d( + const float *in, + float *out, + float *workspace, + std::ptrdiff_t nz, + std::ptrdiff_t ny, + std::ptrdiff_t nx, + const Kernel1D &kz, + const Kernel1D &ky, + const Kernel1D &kx +) { + // Axis 0 (Z), strided: (1, nz, ny*nx). in -> out. + convolve_axis_strided(in, out, 1, nz, ny * nx, kz); + // Axis 1 (Y), strided: (nz, ny, nx). out -> workspace. + convolve_axis_strided(out, workspace, nz, ny, nx, ky); + // Axis 2 (X), contiguous: (nz*ny, nx). workspace -> out. + convolve_axis_x(workspace, out, nz * ny, nx, kx); +} + +// --------------------------------------------------------------------------- +// Public composite filters. All operate on float32 C-contiguous buffers in +// NumPy axis order: 2D = (ny, nx), 3D = (nz, ny, nx). `sigma_*` and `order_*` +// are per-axis; `window_ratio = 0` selects the default kernel radius +// (ceil((3 + 0.5 * order) * sigma) per axis). +// --------------------------------------------------------------------------- + +inline void gaussian_smoothing_2d( + const float *in, + float *out, + std::ptrdiff_t ny, + std::ptrdiff_t nx, + double sigma_y, + double sigma_x, + double window_ratio +) { + const auto ky = gaussian_kernel(sigma_y, 0, window_ratio); + const auto kx = gaussian_kernel(sigma_x, 0, window_ratio); + std::vector workspace(static_cast(ny * nx)); + gaussian_separable_2d(in, out, workspace.data(), ny, nx, ky, kx); +} + +inline void gaussian_smoothing_3d( + const float *in, + float *out, + std::ptrdiff_t nz, + std::ptrdiff_t ny, + std::ptrdiff_t nx, + double sigma_z, + double sigma_y, + double sigma_x, + double window_ratio +) { + const auto kz = gaussian_kernel(sigma_z, 0, window_ratio); + const auto ky = gaussian_kernel(sigma_y, 0, window_ratio); + const auto kx = gaussian_kernel(sigma_x, 0, window_ratio); + std::vector workspace(static_cast(nz * ny * nx)); + gaussian_separable_3d(in, out, workspace.data(), nz, ny, nx, kz, ky, kx); +} + +inline void gaussian_derivative_2d( + const float *in, + float *out, + std::ptrdiff_t ny, + std::ptrdiff_t nx, + double sigma_y, + double sigma_x, + int order_y, + int order_x, + double window_ratio +) { + const auto ky = gaussian_kernel(sigma_y, order_y, window_ratio); + const auto kx = gaussian_kernel(sigma_x, order_x, window_ratio); + std::vector workspace(static_cast(ny * nx)); + gaussian_separable_2d(in, out, workspace.data(), ny, nx, ky, kx); +} + +inline void gaussian_derivative_3d( + const float *in, + float *out, + std::ptrdiff_t nz, + std::ptrdiff_t ny, + std::ptrdiff_t nx, + double sigma_z, + double sigma_y, + double sigma_x, + int order_z, + int order_y, + int order_x, + double window_ratio +) { + const auto kz = gaussian_kernel(sigma_z, order_z, window_ratio); + const auto ky = gaussian_kernel(sigma_y, order_y, window_ratio); + const auto kx = gaussian_kernel(sigma_x, order_x, window_ratio); + std::vector workspace(static_cast(nz * ny * nx)); + gaussian_separable_3d(in, out, workspace.data(), nz, ny, nx, kz, ky, kx); +} + +inline void gaussian_gradient_magnitude_2d( + const float *in, + float *out, + std::ptrdiff_t ny, + std::ptrdiff_t nx, + double sigma_y, + double sigma_x, + double window_ratio +) { + const std::ptrdiff_t n = ny * nx; + const auto ky0 = gaussian_kernel(sigma_y, 0, window_ratio); + const auto kx0 = gaussian_kernel(sigma_x, 0, window_ratio); + const auto ky1 = gaussian_kernel(sigma_y, 1, window_ratio); + const auto kx1 = gaussian_kernel(sigma_x, 1, window_ratio); + + std::vector work1(static_cast(n)); + std::vector work2(static_cast(n)); + + // d/dy + gaussian_separable_2d(in, work1.data(), work2.data(), ny, nx, ky1, kx0); + for (std::ptrdiff_t i = 0; i < n; ++i) out[i] = work1[i] * work1[i]; + + // d/dx + gaussian_separable_2d(in, work1.data(), work2.data(), ny, nx, ky0, kx1); + for (std::ptrdiff_t i = 0; i < n; ++i) out[i] += work1[i] * work1[i]; + + for (std::ptrdiff_t i = 0; i < n; ++i) out[i] = std::sqrt(out[i]); +} + +inline void gaussian_gradient_magnitude_3d( + const float *in, + float *out, + std::ptrdiff_t nz, + std::ptrdiff_t ny, + std::ptrdiff_t nx, + double sigma_z, + double sigma_y, + double sigma_x, + double window_ratio +) { + const std::ptrdiff_t n = nz * ny * nx; + const auto kz0 = gaussian_kernel(sigma_z, 0, window_ratio); + const auto ky0 = gaussian_kernel(sigma_y, 0, window_ratio); + const auto kx0 = gaussian_kernel(sigma_x, 0, window_ratio); + const auto kz1 = gaussian_kernel(sigma_z, 1, window_ratio); + const auto ky1 = gaussian_kernel(sigma_y, 1, window_ratio); + const auto kx1 = gaussian_kernel(sigma_x, 1, window_ratio); + + std::vector work1(static_cast(n)); + std::vector work2(static_cast(n)); + + gaussian_separable_3d(in, work1.data(), work2.data(), nz, ny, nx, kz1, ky0, kx0); + for (std::ptrdiff_t i = 0; i < n; ++i) out[i] = work1[i] * work1[i]; + + gaussian_separable_3d(in, work1.data(), work2.data(), nz, ny, nx, kz0, ky1, kx0); + for (std::ptrdiff_t i = 0; i < n; ++i) out[i] += work1[i] * work1[i]; + + gaussian_separable_3d(in, work1.data(), work2.data(), nz, ny, nx, kz0, ky0, kx1); + for (std::ptrdiff_t i = 0; i < n; ++i) out[i] += work1[i] * work1[i]; + + for (std::ptrdiff_t i = 0; i < n; ++i) out[i] = std::sqrt(out[i]); +} + +inline void laplacian_of_gaussian_2d( + const float *in, + float *out, + std::ptrdiff_t ny, + std::ptrdiff_t nx, + double sigma_y, + double sigma_x, + double window_ratio +) { + const std::ptrdiff_t n = ny * nx; + const auto ky0 = gaussian_kernel(sigma_y, 0, window_ratio); + const auto kx0 = gaussian_kernel(sigma_x, 0, window_ratio); + const auto ky2 = gaussian_kernel(sigma_y, 2, window_ratio); + const auto kx2 = gaussian_kernel(sigma_x, 2, window_ratio); + + std::vector work1(static_cast(n)); + std::vector work2(static_cast(n)); + + // d²/dy² + gaussian_separable_2d(in, work1.data(), work2.data(), ny, nx, ky2, kx0); + for (std::ptrdiff_t i = 0; i < n; ++i) out[i] = work1[i]; + + // d²/dx² + gaussian_separable_2d(in, work1.data(), work2.data(), ny, nx, ky0, kx2); + for (std::ptrdiff_t i = 0; i < n; ++i) out[i] += work1[i]; +} + +inline void laplacian_of_gaussian_3d( + const float *in, + float *out, + std::ptrdiff_t nz, + std::ptrdiff_t ny, + std::ptrdiff_t nx, + double sigma_z, + double sigma_y, + double sigma_x, + double window_ratio +) { + const std::ptrdiff_t n = nz * ny * nx; + const auto kz0 = gaussian_kernel(sigma_z, 0, window_ratio); + const auto ky0 = gaussian_kernel(sigma_y, 0, window_ratio); + const auto kx0 = gaussian_kernel(sigma_x, 0, window_ratio); + const auto kz2 = gaussian_kernel(sigma_z, 2, window_ratio); + const auto ky2 = gaussian_kernel(sigma_y, 2, window_ratio); + const auto kx2 = gaussian_kernel(sigma_x, 2, window_ratio); + + std::vector work1(static_cast(n)); + std::vector work2(static_cast(n)); + + gaussian_separable_3d(in, work1.data(), work2.data(), nz, ny, nx, kz2, ky0, kx0); + for (std::ptrdiff_t i = 0; i < n; ++i) out[i] = work1[i]; + + gaussian_separable_3d(in, work1.data(), work2.data(), nz, ny, nx, kz0, ky2, kx0); + for (std::ptrdiff_t i = 0; i < n; ++i) out[i] += work1[i]; + + gaussian_separable_3d(in, work1.data(), work2.data(), nz, ny, nx, kz0, ky0, kx2); + for (std::ptrdiff_t i = 0; i < n; ++i) out[i] += work1[i]; +} + +// --------------------------------------------------------------------------- +// Hessian-of-Gaussian eigenvalues. `out` has trailing-axis layout: in 2D, the +// output buffer holds ny*nx*2 floats laid out so out[2*i + 0] = largest +// eigenvalue, out[2*i + 1] = smallest. In 3D similarly with stride 3. +// --------------------------------------------------------------------------- + +inline void hessian_of_gaussian_eigenvalues_2d( + const float *in, + float *out, + std::ptrdiff_t ny, + std::ptrdiff_t nx, + double sigma_y, + double sigma_x, + double window_ratio +) { + const std::ptrdiff_t n = ny * nx; + const auto ky0 = gaussian_kernel(sigma_y, 0, window_ratio); + const auto kx0 = gaussian_kernel(sigma_x, 0, window_ratio); + const auto ky1 = gaussian_kernel(sigma_y, 1, window_ratio); + const auto kx1 = gaussian_kernel(sigma_x, 1, window_ratio); + const auto ky2 = gaussian_kernel(sigma_y, 2, window_ratio); + const auto kx2 = gaussian_kernel(sigma_x, 2, window_ratio); + + std::vector work(static_cast(n)); + std::vector hyy(static_cast(n)); + std::vector hyx(static_cast(n)); + std::vector hxx(static_cast(n)); + + // d²/dy² + gaussian_separable_2d(in, hyy.data(), work.data(), ny, nx, ky2, kx0); + // d²/(dy dx) + gaussian_separable_2d(in, hyx.data(), work.data(), ny, nx, ky1, kx1); + // d²/dx² + gaussian_separable_2d(in, hxx.data(), work.data(), ny, nx, ky0, kx2); + + // ev2 sorted descending into interleaved output. + for (std::ptrdiff_t i = 0; i < n; ++i) { + const float a = hyy[i]; + const float b = hyx[i]; + const float c = hxx[i]; + const float half_tr = 0.5f * (a + c); + const float half_diff = 0.5f * (a - c); + const float disc = std::sqrt(half_diff * half_diff + b * b); + out[2 * i + 0] = half_tr + disc; + out[2 * i + 1] = half_tr - disc; + } +} + +inline void hessian_of_gaussian_eigenvalues_3d( + const float *in, + float *out, + std::ptrdiff_t nz, + std::ptrdiff_t ny, + std::ptrdiff_t nx, + double sigma_z, + double sigma_y, + double sigma_x, + double window_ratio +) { + const std::ptrdiff_t n = nz * ny * nx; + const auto kz0 = gaussian_kernel(sigma_z, 0, window_ratio); + const auto ky0 = gaussian_kernel(sigma_y, 0, window_ratio); + const auto kx0 = gaussian_kernel(sigma_x, 0, window_ratio); + const auto kz1 = gaussian_kernel(sigma_z, 1, window_ratio); + const auto ky1 = gaussian_kernel(sigma_y, 1, window_ratio); + const auto kx1 = gaussian_kernel(sigma_x, 1, window_ratio); + const auto kz2 = gaussian_kernel(sigma_z, 2, window_ratio); + const auto ky2 = gaussian_kernel(sigma_y, 2, window_ratio); + const auto kx2 = gaussian_kernel(sigma_x, 2, window_ratio); + + std::vector work(static_cast(n)); + std::vector hzz(static_cast(n)); + std::vector hzy(static_cast(n)); + std::vector hzx(static_cast(n)); + std::vector hyy(static_cast(n)); + std::vector hyx(static_cast(n)); + std::vector hxx(static_cast(n)); + + gaussian_separable_3d(in, hzz.data(), work.data(), nz, ny, nx, kz2, ky0, kx0); + gaussian_separable_3d(in, hzy.data(), work.data(), nz, ny, nx, kz1, ky1, kx0); + gaussian_separable_3d(in, hzx.data(), work.data(), nz, ny, nx, kz1, ky0, kx1); + gaussian_separable_3d(in, hyy.data(), work.data(), nz, ny, nx, kz0, ky2, kx0); + gaussian_separable_3d(in, hyx.data(), work.data(), nz, ny, nx, kz0, ky1, kx1); + gaussian_separable_3d(in, hxx.data(), work.data(), nz, ny, nx, kz0, ky0, kx2); + + // ev3 sorted descending into interleaved output. + for (std::ptrdiff_t i = 0; i < n; ++i) { + float e0; + float e1; + float e2; + detail::ev3_one_descending( + hzz[i], hzy[i], hzx[i], hyy[i], hyx[i], hxx[i], + e0, e1, e2 + ); + out[3 * i + 0] = e0; + out[3 * i + 1] = e1; + out[3 * i + 2] = e2; + } +} + +// --------------------------------------------------------------------------- +// Structure-tensor eigenvalues. Two-scale: first take first-order Gaussian +// derivatives at sigma_inner, form the outer products, smooth them with +// sigma_outer, then compute eigenvalues of the resulting symmetric tensor. +// Output layout matches the Hessian variants (trailing axis size N). +// --------------------------------------------------------------------------- + +inline void structure_tensor_eigenvalues_2d( + const float *in, + float *out, + std::ptrdiff_t ny, + std::ptrdiff_t nx, + double sigma_inner_y, + double sigma_inner_x, + double sigma_outer_y, + double sigma_outer_x, + double window_ratio +) { + const std::ptrdiff_t n = ny * nx; + + const auto kiy0 = gaussian_kernel(sigma_inner_y, 0, window_ratio); + const auto kix0 = gaussian_kernel(sigma_inner_x, 0, window_ratio); + const auto kiy1 = gaussian_kernel(sigma_inner_y, 1, window_ratio); + const auto kix1 = gaussian_kernel(sigma_inner_x, 1, window_ratio); + const auto koy0 = gaussian_kernel(sigma_outer_y, 0, window_ratio); + const auto kox0 = gaussian_kernel(sigma_outer_x, 0, window_ratio); + + std::vector work(static_cast(n)); + std::vector gy(static_cast(n)); + std::vector gx(static_cast(n)); + std::vector tmp(static_cast(n)); + std::vector syy(static_cast(n)); + std::vector syx(static_cast(n)); + std::vector sxx(static_cast(n)); + + // First-order partials at sigma_inner. + gaussian_separable_2d(in, gy.data(), work.data(), ny, nx, kiy1, kix0); + gaussian_separable_2d(in, gx.data(), work.data(), ny, nx, kiy0, kix1); + + // Outer products, smoothed with sigma_outer. + for (std::ptrdiff_t i = 0; i < n; ++i) tmp[i] = gy[i] * gy[i]; + gaussian_separable_2d(tmp.data(), syy.data(), work.data(), ny, nx, koy0, kox0); + + for (std::ptrdiff_t i = 0; i < n; ++i) tmp[i] = gy[i] * gx[i]; + gaussian_separable_2d(tmp.data(), syx.data(), work.data(), ny, nx, koy0, kox0); + + for (std::ptrdiff_t i = 0; i < n; ++i) tmp[i] = gx[i] * gx[i]; + gaussian_separable_2d(tmp.data(), sxx.data(), work.data(), ny, nx, koy0, kox0); + + for (std::ptrdiff_t i = 0; i < n; ++i) { + const float a = syy[i]; + const float b = syx[i]; + const float c = sxx[i]; + const float half_tr = 0.5f * (a + c); + const float half_diff = 0.5f * (a - c); + const float disc = std::sqrt(half_diff * half_diff + b * b); + out[2 * i + 0] = half_tr + disc; + out[2 * i + 1] = half_tr - disc; + } +} + +inline void structure_tensor_eigenvalues_3d( + const float *in, + float *out, + std::ptrdiff_t nz, + std::ptrdiff_t ny, + std::ptrdiff_t nx, + double sigma_inner_z, + double sigma_inner_y, + double sigma_inner_x, + double sigma_outer_z, + double sigma_outer_y, + double sigma_outer_x, + double window_ratio +) { + const std::ptrdiff_t n = nz * ny * nx; + + const auto kiz0 = gaussian_kernel(sigma_inner_z, 0, window_ratio); + const auto kiy0 = gaussian_kernel(sigma_inner_y, 0, window_ratio); + const auto kix0 = gaussian_kernel(sigma_inner_x, 0, window_ratio); + const auto kiz1 = gaussian_kernel(sigma_inner_z, 1, window_ratio); + const auto kiy1 = gaussian_kernel(sigma_inner_y, 1, window_ratio); + const auto kix1 = gaussian_kernel(sigma_inner_x, 1, window_ratio); + const auto koz0 = gaussian_kernel(sigma_outer_z, 0, window_ratio); + const auto koy0 = gaussian_kernel(sigma_outer_y, 0, window_ratio); + const auto kox0 = gaussian_kernel(sigma_outer_x, 0, window_ratio); + + std::vector work(static_cast(n)); + std::vector gz(static_cast(n)); + std::vector gy(static_cast(n)); + std::vector gx(static_cast(n)); + std::vector tmp(static_cast(n)); + std::vector szz(static_cast(n)); + std::vector szy(static_cast(n)); + std::vector szx(static_cast(n)); + std::vector syy(static_cast(n)); + std::vector syx(static_cast(n)); + std::vector sxx(static_cast(n)); + + gaussian_separable_3d(in, gz.data(), work.data(), nz, ny, nx, kiz1, kiy0, kix0); + gaussian_separable_3d(in, gy.data(), work.data(), nz, ny, nx, kiz0, kiy1, kix0); + gaussian_separable_3d(in, gx.data(), work.data(), nz, ny, nx, kiz0, kiy0, kix1); + + for (std::ptrdiff_t i = 0; i < n; ++i) tmp[i] = gz[i] * gz[i]; + gaussian_separable_3d(tmp.data(), szz.data(), work.data(), nz, ny, nx, koz0, koy0, kox0); + + for (std::ptrdiff_t i = 0; i < n; ++i) tmp[i] = gz[i] * gy[i]; + gaussian_separable_3d(tmp.data(), szy.data(), work.data(), nz, ny, nx, koz0, koy0, kox0); + + for (std::ptrdiff_t i = 0; i < n; ++i) tmp[i] = gz[i] * gx[i]; + gaussian_separable_3d(tmp.data(), szx.data(), work.data(), nz, ny, nx, koz0, koy0, kox0); + + for (std::ptrdiff_t i = 0; i < n; ++i) tmp[i] = gy[i] * gy[i]; + gaussian_separable_3d(tmp.data(), syy.data(), work.data(), nz, ny, nx, koz0, koy0, kox0); + + for (std::ptrdiff_t i = 0; i < n; ++i) tmp[i] = gy[i] * gx[i]; + gaussian_separable_3d(tmp.data(), syx.data(), work.data(), nz, ny, nx, koz0, koy0, kox0); + + for (std::ptrdiff_t i = 0; i < n; ++i) tmp[i] = gx[i] * gx[i]; + gaussian_separable_3d(tmp.data(), sxx.data(), work.data(), nz, ny, nx, koz0, koy0, kox0); + + for (std::ptrdiff_t i = 0; i < n; ++i) { + float e0; + float e1; + float e2; + detail::ev3_one_descending( + szz[i], szy[i], szx[i], syy[i], syx[i], sxx[i], + e0, e1, e2 + ); + out[3 * i + 0] = e0; + out[3 * i + 1] = e1; + out[3 * i + 2] = e2; + } +} + +} // namespace bioimage_cpp::filters diff --git a/include/bioimage_cpp/filters/kernel.hxx b/include/bioimage_cpp/filters/kernel.hxx new file mode 100644 index 0000000..7d1d9ec --- /dev/null +++ b/include/bioimage_cpp/filters/kernel.hxx @@ -0,0 +1,137 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +namespace bioimage_cpp::filters { + +// Half-Gaussian-derivative 1D kernel. +// +// Storage convention: +// half_coefs has length `radius + 1`. half_coefs[0] is the centre tap. +// For symmetric kernels (order 0/2), kernel[k] = kernel[-k] = half_coefs[|k|]. +// For antisymmetric kernels (order 1), kernel[k] = sign(k) * half_coefs[|k|], +// with kernel[0] = 0 (so half_coefs[0] is always 0 for antisymmetric). +// +// Convolution convention (cross-correlation): +// output[i] = sum_{k=-R..R} kernel[k] * input[i + k] +// For symmetric: output[i] = h[0]*x[i] + sum_{k>=1} h[k]*(x[i+k] + x[i-k]) +// For antisymmetric: output[i] = sum_{k>=1} h[k]*(x[i+k] - x[i-k]) +// +// Normalisation targets are chosen so that convolving with the polynomial +// f(x) = x^order / order! yields a constant 1 at the centre, i.e. these +// kernels evaluate the derivative directly (not the derivative scaled by +// sigma^order, in contrast to some other conventions). +struct Kernel1D { + std::vector half_coefs; + int radius = 0; + int order = 0; + bool is_symmetric = true; +}; + +inline int gaussian_kernel_radius(double sigma, int order, double window_ratio) { + const double effective_ratio = + (window_ratio > 0.0) ? window_ratio : (3.0 + 0.5 * static_cast(order)); + const int r = static_cast(std::ceil(effective_ratio * sigma)); + return std::max(1, r); +} + +// Build a 1D Gaussian (derivative) kernel. +// sigma > 0; order in {0, 1, 2}; window_ratio = 0 selects the default radius +// (ceil((3 + 0.5*order)*sigma)). +inline Kernel1D gaussian_kernel(double sigma, int order, double window_ratio = 0.0) { + if (!(sigma > 0.0)) { + throw std::invalid_argument( + "sigma must be positive, got sigma=" + std::to_string(sigma) + ); + } + if (order < 0 || order > 2) { + throw std::invalid_argument( + "order must be 0, 1 or 2, got order=" + std::to_string(order) + ); + } + + const int radius = gaussian_kernel_radius(sigma, order, window_ratio); + const double inv_sigma2 = 1.0 / (sigma * sigma); + const double inv_2sigma2 = 0.5 * inv_sigma2; + + Kernel1D kernel; + kernel.radius = radius; + kernel.order = order; + kernel.is_symmetric = (order != 1); + kernel.half_coefs.assign(static_cast(radius + 1), 0.0f); + + // Generate raw (unnormalised) half-kernel values in double precision. + std::vector raw(static_cast(radius + 1)); + for (int i = 0; i <= radius; ++i) { + const double x = static_cast(i); + const double g = std::exp(-x * x * inv_2sigma2); + switch (order) { + case 0: + raw[i] = g; + break; + case 1: + // Store +x/sigma^2 * g(x) for i>=0 so the antisymmetric loop + // combines as (x[i+k] - x[i-k]) and produces the standard + // positive d/dx convention. + raw[i] = (x * inv_sigma2) * g; + break; + case 2: + raw[i] = (x * x * inv_sigma2 * inv_sigma2 - inv_sigma2) * g; + break; + } + } + + if (order == 0) { + // Sum to 1. + double sum = raw[0]; + for (int i = 1; i <= radius; ++i) sum += 2.0 * raw[i]; + const double scale = 1.0 / sum; + for (int i = 0; i <= radius; ++i) { + kernel.half_coefs[static_cast(i)] = static_cast(raw[i] * scale); + } + } else if (order == 1) { + // Antisymmetric. Normalise so sum_{k>=1} k * h[k] = 0.5, i.e. the + // discrete first derivative of f(x)=x evaluates to 1. + raw[0] = 0.0; // ensure centre tap is exactly zero + double moment = 0.0; + for (int i = 1; i <= radius; ++i) { + moment += static_cast(i) * raw[i]; + } + if (moment == 0.0) { + throw std::runtime_error("degenerate Gaussian first-derivative kernel"); + } + const double scale = 0.5 / moment; + for (int i = 0; i <= radius; ++i) { + kernel.half_coefs[static_cast(i)] = static_cast(raw[i] * scale); + } + } else { + // order == 2. DC-correct (kernel sums to 0), then normalise so + // sum_{k>=1} k^2 * h[k] = 1, i.e. discrete second derivative of + // f(x) = x^2/2 evaluates to 1. + double sum = raw[0]; + for (int i = 1; i <= radius; ++i) sum += 2.0 * raw[i]; + const double dc = sum / static_cast(2 * radius + 1); + for (int i = 0; i <= radius; ++i) raw[i] -= dc; + + double moment = 0.0; + for (int i = 1; i <= radius; ++i) { + moment += static_cast(i) * static_cast(i) * raw[i]; + } + if (moment == 0.0) { + throw std::runtime_error("degenerate Gaussian second-derivative kernel"); + } + const double scale = 1.0 / moment; + for (int i = 0; i <= radius; ++i) { + kernel.half_coefs[static_cast(i)] = static_cast(raw[i] * scale); + } + } + + return kernel; +} + +} // namespace bioimage_cpp::filters diff --git a/src/bindings/filters.cxx b/src/bindings/filters.cxx new file mode 100644 index 0000000..396cb63 --- /dev/null +++ b/src/bindings/filters.cxx @@ -0,0 +1,565 @@ +#include "filters.hxx" + +#include "bioimage_cpp/filters/gaussian.hxx" + +#include + +#include +#include +#include + +namespace nb = nanobind; + +namespace bioimage_cpp::bindings { +namespace { + +using ConstImage = nb::ndarray; +using Image = nb::ndarray; + +void require_ndim(const ConstImage &image, int expected, const char *function) { + if (static_cast(image.ndim()) != expected) { + throw std::invalid_argument( + std::string(function) + ": image must have ndim=" + std::to_string(expected) + + ", got ndim=" + std::to_string(image.ndim()) + ); + } +} + +void require_positive_sigma(double sigma, const char *name, const char *function) { + if (!(sigma > 0.0)) { + throw std::invalid_argument( + std::string(function) + ": " + name + " must be positive, got " + + std::to_string(sigma) + ); + } +} + +void require_order(int order, const char *name, const char *function) { + if (order < 0 || order > 2) { + throw std::invalid_argument( + std::string(function) + ": " + name + " must be 0, 1 or 2, got " + + std::to_string(order) + ); + } +} + +void require_non_negative_window(double window_ratio, const char *function) { + if (window_ratio < 0.0) { + throw std::invalid_argument( + std::string(function) + ": window_ratio must be >= 0 (0 selects the " + "default), got " + std::to_string(window_ratio) + ); + } +} + +Image allocate_image(const std::size_t *shape, std::size_t ndim) { + std::size_t total = 1; + for (std::size_t i = 0; i < ndim; ++i) total *= shape[i]; + auto *data = new float[total](); + nb::capsule owner(data, [](void *p) noexcept { delete[] static_cast(p); }); + return Image(data, ndim, shape, owner); +} + +// --------------------------------------------------------------------------- +// gaussian_smoothing +// --------------------------------------------------------------------------- + +Image gaussian_smoothing_2d( + ConstImage image, double sigma_y, double sigma_x, double window_ratio +) { + const char *fn = "gaussian_smoothing_2d"; + require_ndim(image, 2, fn); + require_positive_sigma(sigma_y, "sigma_y", fn); + require_positive_sigma(sigma_x, "sigma_x", fn); + require_non_negative_window(window_ratio, fn); + + const std::size_t ny = image.shape(0); + const std::size_t nx = image.shape(1); + const std::size_t shape[2] = {ny, nx}; + Image out = allocate_image(shape, 2); + + const float *in_ptr = image.data(); + float *out_ptr = out.data(); + { + nb::gil_scoped_release release; + filters::gaussian_smoothing_2d( + in_ptr, out_ptr, + static_cast(ny), + static_cast(nx), + sigma_y, sigma_x, window_ratio + ); + } + return out; +} + +Image gaussian_smoothing_3d( + ConstImage image, + double sigma_z, double sigma_y, double sigma_x, double window_ratio +) { + const char *fn = "gaussian_smoothing_3d"; + require_ndim(image, 3, fn); + require_positive_sigma(sigma_z, "sigma_z", fn); + require_positive_sigma(sigma_y, "sigma_y", fn); + require_positive_sigma(sigma_x, "sigma_x", fn); + require_non_negative_window(window_ratio, fn); + + const std::size_t nz = image.shape(0); + const std::size_t ny = image.shape(1); + const std::size_t nx = image.shape(2); + const std::size_t shape[3] = {nz, ny, nx}; + Image out = allocate_image(shape, 3); + + const float *in_ptr = image.data(); + float *out_ptr = out.data(); + { + nb::gil_scoped_release release; + filters::gaussian_smoothing_3d( + in_ptr, out_ptr, + static_cast(nz), + static_cast(ny), + static_cast(nx), + sigma_z, sigma_y, sigma_x, window_ratio + ); + } + return out; +} + +// --------------------------------------------------------------------------- +// gaussian_derivative +// --------------------------------------------------------------------------- + +Image gaussian_derivative_2d( + ConstImage image, + double sigma_y, double sigma_x, + int order_y, int order_x, + double window_ratio +) { + const char *fn = "gaussian_derivative_2d"; + require_ndim(image, 2, fn); + require_positive_sigma(sigma_y, "sigma_y", fn); + require_positive_sigma(sigma_x, "sigma_x", fn); + require_order(order_y, "order_y", fn); + require_order(order_x, "order_x", fn); + require_non_negative_window(window_ratio, fn); + + const std::size_t ny = image.shape(0); + const std::size_t nx = image.shape(1); + const std::size_t shape[2] = {ny, nx}; + Image out = allocate_image(shape, 2); + + const float *in_ptr = image.data(); + float *out_ptr = out.data(); + { + nb::gil_scoped_release release; + filters::gaussian_derivative_2d( + in_ptr, out_ptr, + static_cast(ny), + static_cast(nx), + sigma_y, sigma_x, order_y, order_x, window_ratio + ); + } + return out; +} + +Image gaussian_derivative_3d( + ConstImage image, + double sigma_z, double sigma_y, double sigma_x, + int order_z, int order_y, int order_x, + double window_ratio +) { + const char *fn = "gaussian_derivative_3d"; + require_ndim(image, 3, fn); + require_positive_sigma(sigma_z, "sigma_z", fn); + require_positive_sigma(sigma_y, "sigma_y", fn); + require_positive_sigma(sigma_x, "sigma_x", fn); + require_order(order_z, "order_z", fn); + require_order(order_y, "order_y", fn); + require_order(order_x, "order_x", fn); + require_non_negative_window(window_ratio, fn); + + const std::size_t nz = image.shape(0); + const std::size_t ny = image.shape(1); + const std::size_t nx = image.shape(2); + const std::size_t shape[3] = {nz, ny, nx}; + Image out = allocate_image(shape, 3); + + const float *in_ptr = image.data(); + float *out_ptr = out.data(); + { + nb::gil_scoped_release release; + filters::gaussian_derivative_3d( + in_ptr, out_ptr, + static_cast(nz), + static_cast(ny), + static_cast(nx), + sigma_z, sigma_y, sigma_x, + order_z, order_y, order_x, + window_ratio + ); + } + return out; +} + +// --------------------------------------------------------------------------- +// gradient magnitude +// --------------------------------------------------------------------------- + +Image gaussian_gradient_magnitude_2d( + ConstImage image, double sigma_y, double sigma_x, double window_ratio +) { + const char *fn = "gaussian_gradient_magnitude_2d"; + require_ndim(image, 2, fn); + require_positive_sigma(sigma_y, "sigma_y", fn); + require_positive_sigma(sigma_x, "sigma_x", fn); + require_non_negative_window(window_ratio, fn); + + const std::size_t ny = image.shape(0); + const std::size_t nx = image.shape(1); + const std::size_t shape[2] = {ny, nx}; + Image out = allocate_image(shape, 2); + + const float *in_ptr = image.data(); + float *out_ptr = out.data(); + { + nb::gil_scoped_release release; + filters::gaussian_gradient_magnitude_2d( + in_ptr, out_ptr, + static_cast(ny), + static_cast(nx), + sigma_y, sigma_x, window_ratio + ); + } + return out; +} + +Image gaussian_gradient_magnitude_3d( + ConstImage image, + double sigma_z, double sigma_y, double sigma_x, double window_ratio +) { + const char *fn = "gaussian_gradient_magnitude_3d"; + require_ndim(image, 3, fn); + require_positive_sigma(sigma_z, "sigma_z", fn); + require_positive_sigma(sigma_y, "sigma_y", fn); + require_positive_sigma(sigma_x, "sigma_x", fn); + require_non_negative_window(window_ratio, fn); + + const std::size_t nz = image.shape(0); + const std::size_t ny = image.shape(1); + const std::size_t nx = image.shape(2); + const std::size_t shape[3] = {nz, ny, nx}; + Image out = allocate_image(shape, 3); + + const float *in_ptr = image.data(); + float *out_ptr = out.data(); + { + nb::gil_scoped_release release; + filters::gaussian_gradient_magnitude_3d( + in_ptr, out_ptr, + static_cast(nz), + static_cast(ny), + static_cast(nx), + sigma_z, sigma_y, sigma_x, window_ratio + ); + } + return out; +} + +// --------------------------------------------------------------------------- +// Laplacian of Gaussian +// --------------------------------------------------------------------------- + +Image laplacian_of_gaussian_2d( + ConstImage image, double sigma_y, double sigma_x, double window_ratio +) { + const char *fn = "laplacian_of_gaussian_2d"; + require_ndim(image, 2, fn); + require_positive_sigma(sigma_y, "sigma_y", fn); + require_positive_sigma(sigma_x, "sigma_x", fn); + require_non_negative_window(window_ratio, fn); + + const std::size_t ny = image.shape(0); + const std::size_t nx = image.shape(1); + const std::size_t shape[2] = {ny, nx}; + Image out = allocate_image(shape, 2); + + const float *in_ptr = image.data(); + float *out_ptr = out.data(); + { + nb::gil_scoped_release release; + filters::laplacian_of_gaussian_2d( + in_ptr, out_ptr, + static_cast(ny), + static_cast(nx), + sigma_y, sigma_x, window_ratio + ); + } + return out; +} + +Image laplacian_of_gaussian_3d( + ConstImage image, + double sigma_z, double sigma_y, double sigma_x, double window_ratio +) { + const char *fn = "laplacian_of_gaussian_3d"; + require_ndim(image, 3, fn); + require_positive_sigma(sigma_z, "sigma_z", fn); + require_positive_sigma(sigma_y, "sigma_y", fn); + require_positive_sigma(sigma_x, "sigma_x", fn); + require_non_negative_window(window_ratio, fn); + + const std::size_t nz = image.shape(0); + const std::size_t ny = image.shape(1); + const std::size_t nx = image.shape(2); + const std::size_t shape[3] = {nz, ny, nx}; + Image out = allocate_image(shape, 3); + + const float *in_ptr = image.data(); + float *out_ptr = out.data(); + { + nb::gil_scoped_release release; + filters::laplacian_of_gaussian_3d( + in_ptr, out_ptr, + static_cast(nz), + static_cast(ny), + static_cast(nx), + sigma_z, sigma_y, sigma_x, window_ratio + ); + } + return out; +} + +// --------------------------------------------------------------------------- +// Hessian-of-Gaussian eigenvalues. Output shape: input shape + (N,) trailing. +// --------------------------------------------------------------------------- + +Image hessian_of_gaussian_eigenvalues_2d( + ConstImage image, double sigma_y, double sigma_x, double window_ratio +) { + const char *fn = "hessian_of_gaussian_eigenvalues_2d"; + require_ndim(image, 2, fn); + require_positive_sigma(sigma_y, "sigma_y", fn); + require_positive_sigma(sigma_x, "sigma_x", fn); + require_non_negative_window(window_ratio, fn); + + const std::size_t ny = image.shape(0); + const std::size_t nx = image.shape(1); + const std::size_t shape[3] = {ny, nx, 2}; + Image out = allocate_image(shape, 3); + + const float *in_ptr = image.data(); + float *out_ptr = out.data(); + { + nb::gil_scoped_release release; + filters::hessian_of_gaussian_eigenvalues_2d( + in_ptr, out_ptr, + static_cast(ny), + static_cast(nx), + sigma_y, sigma_x, window_ratio + ); + } + return out; +} + +Image hessian_of_gaussian_eigenvalues_3d( + ConstImage image, + double sigma_z, double sigma_y, double sigma_x, double window_ratio +) { + const char *fn = "hessian_of_gaussian_eigenvalues_3d"; + require_ndim(image, 3, fn); + require_positive_sigma(sigma_z, "sigma_z", fn); + require_positive_sigma(sigma_y, "sigma_y", fn); + require_positive_sigma(sigma_x, "sigma_x", fn); + require_non_negative_window(window_ratio, fn); + + const std::size_t nz = image.shape(0); + const std::size_t ny = image.shape(1); + const std::size_t nx = image.shape(2); + const std::size_t shape[4] = {nz, ny, nx, 3}; + Image out = allocate_image(shape, 4); + + const float *in_ptr = image.data(); + float *out_ptr = out.data(); + { + nb::gil_scoped_release release; + filters::hessian_of_gaussian_eigenvalues_3d( + in_ptr, out_ptr, + static_cast(nz), + static_cast(ny), + static_cast(nx), + sigma_z, sigma_y, sigma_x, window_ratio + ); + } + return out; +} + +// --------------------------------------------------------------------------- +// Structure-tensor eigenvalues. Output shape: input shape + (N,) trailing. +// --------------------------------------------------------------------------- + +Image structure_tensor_eigenvalues_2d( + ConstImage image, + double sigma_inner_y, double sigma_inner_x, + double sigma_outer_y, double sigma_outer_x, + double window_ratio +) { + const char *fn = "structure_tensor_eigenvalues_2d"; + require_ndim(image, 2, fn); + require_positive_sigma(sigma_inner_y, "sigma_inner_y", fn); + require_positive_sigma(sigma_inner_x, "sigma_inner_x", fn); + require_positive_sigma(sigma_outer_y, "sigma_outer_y", fn); + require_positive_sigma(sigma_outer_x, "sigma_outer_x", fn); + require_non_negative_window(window_ratio, fn); + + const std::size_t ny = image.shape(0); + const std::size_t nx = image.shape(1); + const std::size_t shape[3] = {ny, nx, 2}; + Image out = allocate_image(shape, 3); + + const float *in_ptr = image.data(); + float *out_ptr = out.data(); + { + nb::gil_scoped_release release; + filters::structure_tensor_eigenvalues_2d( + in_ptr, out_ptr, + static_cast(ny), + static_cast(nx), + sigma_inner_y, sigma_inner_x, + sigma_outer_y, sigma_outer_x, + window_ratio + ); + } + return out; +} + +Image structure_tensor_eigenvalues_3d( + ConstImage image, + double sigma_inner_z, double sigma_inner_y, double sigma_inner_x, + double sigma_outer_z, double sigma_outer_y, double sigma_outer_x, + double window_ratio +) { + const char *fn = "structure_tensor_eigenvalues_3d"; + require_ndim(image, 3, fn); + require_positive_sigma(sigma_inner_z, "sigma_inner_z", fn); + require_positive_sigma(sigma_inner_y, "sigma_inner_y", fn); + require_positive_sigma(sigma_inner_x, "sigma_inner_x", fn); + require_positive_sigma(sigma_outer_z, "sigma_outer_z", fn); + require_positive_sigma(sigma_outer_y, "sigma_outer_y", fn); + require_positive_sigma(sigma_outer_x, "sigma_outer_x", fn); + require_non_negative_window(window_ratio, fn); + + const std::size_t nz = image.shape(0); + const std::size_t ny = image.shape(1); + const std::size_t nx = image.shape(2); + const std::size_t shape[4] = {nz, ny, nx, 3}; + Image out = allocate_image(shape, 4); + + const float *in_ptr = image.data(); + float *out_ptr = out.data(); + { + nb::gil_scoped_release release; + filters::structure_tensor_eigenvalues_3d( + in_ptr, out_ptr, + static_cast(nz), + static_cast(ny), + static_cast(nx), + sigma_inner_z, sigma_inner_y, sigma_inner_x, + sigma_outer_z, sigma_outer_y, sigma_outer_x, + window_ratio + ); + } + return out; +} + +} // namespace + +void bind_filters(nb::module_ &m) { + m.def( + "_gaussian_smoothing_2d_float32", &gaussian_smoothing_2d, + nb::arg("image"), nb::arg("sigma_y"), nb::arg("sigma_x"), + nb::arg("window_ratio") = 0.0, + "2D Gaussian smoothing on a float32 (ny, nx) image with anisotropic sigma." + ); + m.def( + "_gaussian_smoothing_3d_float32", &gaussian_smoothing_3d, + nb::arg("image"), + nb::arg("sigma_z"), nb::arg("sigma_y"), nb::arg("sigma_x"), + nb::arg("window_ratio") = 0.0, + "3D Gaussian smoothing on a float32 (nz, ny, nx) image with anisotropic sigma." + ); + m.def( + "_gaussian_derivative_2d_float32", &gaussian_derivative_2d, + nb::arg("image"), nb::arg("sigma_y"), nb::arg("sigma_x"), + nb::arg("order_y"), nb::arg("order_x"), + nb::arg("window_ratio") = 0.0, + "2D Gaussian derivative on a float32 (ny, nx) image with per-axis order." + ); + m.def( + "_gaussian_derivative_3d_float32", &gaussian_derivative_3d, + nb::arg("image"), + nb::arg("sigma_z"), nb::arg("sigma_y"), nb::arg("sigma_x"), + nb::arg("order_z"), nb::arg("order_y"), nb::arg("order_x"), + nb::arg("window_ratio") = 0.0, + "3D Gaussian derivative on a float32 (nz, ny, nx) image with per-axis order." + ); + m.def( + "_gaussian_gradient_magnitude_2d_float32", &gaussian_gradient_magnitude_2d, + nb::arg("image"), nb::arg("sigma_y"), nb::arg("sigma_x"), + nb::arg("window_ratio") = 0.0, + "Gradient magnitude of a Gaussian-smoothed 2D float32 image." + ); + m.def( + "_gaussian_gradient_magnitude_3d_float32", &gaussian_gradient_magnitude_3d, + nb::arg("image"), + nb::arg("sigma_z"), nb::arg("sigma_y"), nb::arg("sigma_x"), + nb::arg("window_ratio") = 0.0, + "Gradient magnitude of a Gaussian-smoothed 3D float32 image." + ); + m.def( + "_laplacian_of_gaussian_2d_float32", &laplacian_of_gaussian_2d, + nb::arg("image"), nb::arg("sigma_y"), nb::arg("sigma_x"), + nb::arg("window_ratio") = 0.0, + "Laplacian of Gaussian on a 2D float32 image." + ); + m.def( + "_laplacian_of_gaussian_3d_float32", &laplacian_of_gaussian_3d, + nb::arg("image"), + nb::arg("sigma_z"), nb::arg("sigma_y"), nb::arg("sigma_x"), + nb::arg("window_ratio") = 0.0, + "Laplacian of Gaussian on a 3D float32 image." + ); + m.def( + "_hessian_of_gaussian_eigenvalues_2d_float32", &hessian_of_gaussian_eigenvalues_2d, + nb::arg("image"), nb::arg("sigma_y"), nb::arg("sigma_x"), + nb::arg("window_ratio") = 0.0, + "Eigenvalues of the Hessian of Gaussian on a 2D float32 image. " + "Output shape: (ny, nx, 2), sorted descending along the trailing axis." + ); + m.def( + "_hessian_of_gaussian_eigenvalues_3d_float32", &hessian_of_gaussian_eigenvalues_3d, + nb::arg("image"), + nb::arg("sigma_z"), nb::arg("sigma_y"), nb::arg("sigma_x"), + nb::arg("window_ratio") = 0.0, + "Eigenvalues of the Hessian of Gaussian on a 3D float32 image. " + "Output shape: (nz, ny, nx, 3), sorted descending along the trailing axis." + ); + m.def( + "_structure_tensor_eigenvalues_2d_float32", &structure_tensor_eigenvalues_2d, + nb::arg("image"), + nb::arg("sigma_inner_y"), nb::arg("sigma_inner_x"), + nb::arg("sigma_outer_y"), nb::arg("sigma_outer_x"), + nb::arg("window_ratio") = 0.0, + "Eigenvalues of the structure tensor on a 2D float32 image. " + "Output shape: (ny, nx, 2), sorted descending along the trailing axis." + ); + m.def( + "_structure_tensor_eigenvalues_3d_float32", &structure_tensor_eigenvalues_3d, + nb::arg("image"), + nb::arg("sigma_inner_z"), nb::arg("sigma_inner_y"), nb::arg("sigma_inner_x"), + nb::arg("sigma_outer_z"), nb::arg("sigma_outer_y"), nb::arg("sigma_outer_x"), + nb::arg("window_ratio") = 0.0, + "Eigenvalues of the structure tensor on a 3D float32 image. " + "Output shape: (nz, ny, nx, 3), sorted descending along the trailing axis." + ); +} + +} // namespace bioimage_cpp::bindings diff --git a/src/bindings/filters.hxx b/src/bindings/filters.hxx new file mode 100644 index 0000000..02e442f --- /dev/null +++ b/src/bindings/filters.hxx @@ -0,0 +1,9 @@ +#pragma once + +#include + +namespace bioimage_cpp::bindings { + +void bind_filters(nanobind::module_ &m); + +} // namespace bioimage_cpp::bindings diff --git a/src/bindings/module.cxx b/src/bindings/module.cxx index cee23a9..a8549f5 100644 --- a/src/bindings/module.cxx +++ b/src/bindings/module.cxx @@ -1,4 +1,5 @@ #include "blocking.hxx" +#include "filters.hxx" #include "graph.hxx" #include "ground_truth.hxx" #include "segmentation.hxx" @@ -11,6 +12,7 @@ namespace nb = nanobind; NB_MODULE(_core, m) { m.doc() = "C++ extension module for bioimage_cpp."; bioimage_cpp::bindings::bind_blocking(m); + bioimage_cpp::bindings::bind_filters(m); bioimage_cpp::bindings::bind_graph(m); bioimage_cpp::bindings::bind_ground_truth(m); bioimage_cpp::bindings::bind_segmentation(m); diff --git a/src/bioimage_cpp/__init__.py b/src/bioimage_cpp/__init__.py index f8a5e6f..2636ebc 100644 --- a/src/bioimage_cpp/__init__.py +++ b/src/bioimage_cpp/__init__.py @@ -2,6 +2,7 @@ from ._version import __version__ from ._core import Block, Blocking, BlockWithHalo +from . import filters from . import graph from . import ground_truth from . import segmentation @@ -12,6 +13,7 @@ "Block", "Blocking", "BlockWithHalo", + "filters", "graph", "ground_truth", "segmentation", diff --git a/src/bioimage_cpp/filters/__init__.py b/src/bioimage_cpp/filters/__init__.py new file mode 100644 index 0000000..716d3fe --- /dev/null +++ b/src/bioimage_cpp/filters/__init__.py @@ -0,0 +1,20 @@ +"""Image filters: separable Gaussian-family derivatives, gradient magnitude, +Laplacian of Gaussian, Hessian and structure-tensor eigenvalues.""" + +from ._filters import ( + gaussian_derivative, + gaussian_gradient_magnitude, + gaussian_smoothing, + hessian_of_gaussian_eigenvalues, + laplacian_of_gaussian, + structure_tensor_eigenvalues, +) + +__all__ = [ + "gaussian_smoothing", + "gaussian_derivative", + "gaussian_gradient_magnitude", + "laplacian_of_gaussian", + "hessian_of_gaussian_eigenvalues", + "structure_tensor_eigenvalues", +] diff --git a/src/bioimage_cpp/filters/_filters.py b/src/bioimage_cpp/filters/_filters.py new file mode 100644 index 0000000..1e1b56e --- /dev/null +++ b/src/bioimage_cpp/filters/_filters.py @@ -0,0 +1,236 @@ +"""Python wrappers for the ``_core`` filter bindings. + +Conventions +----------- +* Input is a 2D or 3D NumPy array (single channel). Supported dtypes: + ``float32``, ``float64``, ``uint8``, ``uint16``. Non-``float32`` inputs are + cast to ``float32`` for the kernel; ``float64`` results are cast back. +* ``sigma`` (and ``order``) accept a scalar or a per-axis sequence of length + ``image.ndim``. Sigmas must be positive; orders must be in ``{0, 1, 2}``. +* ``window_size`` controls the kernel radius as + ``radius = ceil(window_size * sigma)`` (per axis), matching the + ``fastfilters`` / ``vigra`` parameter. ``0`` selects the default + ``3 + 0.5 * order``. +* Axis order is NumPy native: ``(ny, nx)`` for 2D, ``(nz, ny, nx)`` for 3D. +* Eigenvalue functions return an array with a trailing axis of size + ``image.ndim``, sorted descending. +""" + +from __future__ import annotations + +from collections.abc import Sequence + +import numpy as np + +from .. import _core + +_FLOAT_INPUT_DTYPES = (np.float32, np.float64) +_INT_INPUT_DTYPES = (np.uint8, np.uint16) + + +def _prepare_input(image: np.ndarray, function: str) -> tuple[np.ndarray, np.dtype]: + """Validate ndim and cast to contiguous float32. Return (float32 view, + desired output dtype).""" + if not isinstance(image, np.ndarray): + image = np.asarray(image) + if image.ndim not in (2, 3): + raise ValueError( + f"{function}: image must be 2D or 3D, got ndim={image.ndim}" + ) + if image.dtype == np.float32: + return np.ascontiguousarray(image), np.dtype(np.float32) + if image.dtype == np.float64: + return np.ascontiguousarray(image, dtype=np.float32), np.dtype(np.float64) + if image.dtype in (np.dtype(t) for t in _INT_INPUT_DTYPES): + return np.ascontiguousarray(image, dtype=np.float32), np.dtype(np.float32) + raise TypeError( + f"{function}: image dtype must be one of (float32, float64, uint8, " + f"uint16), got dtype={image.dtype}" + ) + + +def _broadcast_per_axis( + value: float | Sequence[float], + ndim: int, + name: str, + function: str, +) -> tuple[float, ...]: + if np.isscalar(value): + return (float(value),) * ndim + seq = tuple(float(v) for v in value) + if len(seq) != ndim: + raise ValueError( + f"{function}: {name} must be a scalar or a sequence of length " + f"{ndim}, got length {len(seq)}" + ) + return seq + + +def _broadcast_order( + value: int | Sequence[int], + ndim: int, + function: str, +) -> tuple[int, ...]: + if np.isscalar(value): + return (int(value),) * ndim + seq = tuple(int(v) for v in value) + if len(seq) != ndim: + raise ValueError( + f"{function}: order must be a scalar or a sequence of length " + f"{ndim}, got length {len(seq)}" + ) + return seq + + +def _finalise(result: np.ndarray, out_dtype: np.dtype) -> np.ndarray: + if result.dtype == out_dtype: + return result + return result.astype(out_dtype, copy=False) + + +def gaussian_smoothing( + image: np.ndarray, + sigma: float | Sequence[float], + *, + window_size: float = 0.0, +) -> np.ndarray: + """Gaussian smoothing of a 2D or 3D scalar image.""" + function = "gaussian_smoothing" + prepared, out_dtype = _prepare_input(image, function) + sigmas = _broadcast_per_axis(sigma, prepared.ndim, "sigma", function) + if prepared.ndim == 2: + result = _core._gaussian_smoothing_2d_float32( + prepared, sigmas[0], sigmas[1], float(window_size) + ) + else: + result = _core._gaussian_smoothing_3d_float32( + prepared, sigmas[0], sigmas[1], sigmas[2], float(window_size) + ) + return _finalise(result, out_dtype) + + +def gaussian_derivative( + image: np.ndarray, + sigma: float | Sequence[float], + order: int | Sequence[int], + *, + window_size: float = 0.0, +) -> np.ndarray: + """Gaussian derivative of a 2D or 3D scalar image with per-axis order.""" + function = "gaussian_derivative" + prepared, out_dtype = _prepare_input(image, function) + sigmas = _broadcast_per_axis(sigma, prepared.ndim, "sigma", function) + orders = _broadcast_order(order, prepared.ndim, function) + if prepared.ndim == 2: + result = _core._gaussian_derivative_2d_float32( + prepared, sigmas[0], sigmas[1], + orders[0], orders[1], float(window_size), + ) + else: + result = _core._gaussian_derivative_3d_float32( + prepared, sigmas[0], sigmas[1], sigmas[2], + orders[0], orders[1], orders[2], float(window_size), + ) + return _finalise(result, out_dtype) + + +def gaussian_gradient_magnitude( + image: np.ndarray, + sigma: float | Sequence[float], + *, + window_size: float = 0.0, +) -> np.ndarray: + """L2 norm of the Gaussian gradient of a 2D or 3D scalar image.""" + function = "gaussian_gradient_magnitude" + prepared, out_dtype = _prepare_input(image, function) + sigmas = _broadcast_per_axis(sigma, prepared.ndim, "sigma", function) + if prepared.ndim == 2: + result = _core._gaussian_gradient_magnitude_2d_float32( + prepared, sigmas[0], sigmas[1], float(window_size) + ) + else: + result = _core._gaussian_gradient_magnitude_3d_float32( + prepared, sigmas[0], sigmas[1], sigmas[2], float(window_size) + ) + return _finalise(result, out_dtype) + + +def laplacian_of_gaussian( + image: np.ndarray, + sigma: float | Sequence[float], + *, + window_size: float = 0.0, +) -> np.ndarray: + """Laplacian of Gaussian (sum of second derivatives) of a 2D or 3D scalar + image.""" + function = "laplacian_of_gaussian" + prepared, out_dtype = _prepare_input(image, function) + sigmas = _broadcast_per_axis(sigma, prepared.ndim, "sigma", function) + if prepared.ndim == 2: + result = _core._laplacian_of_gaussian_2d_float32( + prepared, sigmas[0], sigmas[1], float(window_size) + ) + else: + result = _core._laplacian_of_gaussian_3d_float32( + prepared, sigmas[0], sigmas[1], sigmas[2], float(window_size) + ) + return _finalise(result, out_dtype) + + +def hessian_of_gaussian_eigenvalues( + image: np.ndarray, + sigma: float | Sequence[float], + *, + window_size: float = 0.0, +) -> np.ndarray: + """Eigenvalues of the Hessian of Gaussian. + + Output shape: ``image.shape + (image.ndim,)``, sorted descending along the + trailing axis (largest absolute curvature first only when all eigenvalues + have the same sign; otherwise simply largest signed value first). + """ + function = "hessian_of_gaussian_eigenvalues" + prepared, out_dtype = _prepare_input(image, function) + sigmas = _broadcast_per_axis(sigma, prepared.ndim, "sigma", function) + if prepared.ndim == 2: + result = _core._hessian_of_gaussian_eigenvalues_2d_float32( + prepared, sigmas[0], sigmas[1], float(window_size) + ) + else: + result = _core._hessian_of_gaussian_eigenvalues_3d_float32( + prepared, sigmas[0], sigmas[1], sigmas[2], float(window_size) + ) + return _finalise(result, out_dtype) + + +def structure_tensor_eigenvalues( + image: np.ndarray, + inner_sigma: float | Sequence[float], + outer_sigma: float | Sequence[float], + *, + window_size: float = 0.0, +) -> np.ndarray: + """Eigenvalues of the structure tensor. + + Output shape: ``image.shape + (image.ndim,)``, sorted descending along the + trailing axis. + """ + function = "structure_tensor_eigenvalues" + prepared, out_dtype = _prepare_input(image, function) + inner = _broadcast_per_axis(inner_sigma, prepared.ndim, "inner_sigma", function) + outer = _broadcast_per_axis(outer_sigma, prepared.ndim, "outer_sigma", function) + if prepared.ndim == 2: + result = _core._structure_tensor_eigenvalues_2d_float32( + prepared, + inner[0], inner[1], + outer[0], outer[1], + float(window_size), + ) + else: + result = _core._structure_tensor_eigenvalues_3d_float32( + prepared, + inner[0], inner[1], inner[2], + outer[0], outer[1], outer[2], + float(window_size), + ) + return _finalise(result, out_dtype) diff --git a/tests/test_filters.py b/tests/test_filters.py new file mode 100644 index 0000000..13d6dc8 --- /dev/null +++ b/tests/test_filters.py @@ -0,0 +1,271 @@ +"""Tests for ``bioimage_cpp.filters``. + +The C++ kernels are validated against ``scipy.ndimage`` reference filters with +``mode="mirror"`` (matching our boundary handling). float32 tolerance is +``atol=1e-3`` for composite filters, looser for eigenvalues because of the +trigonometric closed-form rounding. +""" + +import numpy as np +import pytest +from scipy import ndimage + +import bioimage_cpp.filters as bf + + +SHAPES_2D = [(7, 11), (32, 32), (48, 64)] +SHAPES_3D = [(5, 7, 11), (8, 12, 16)] +SIGMAS = [0.7, 1.5, 3.0] + + +def _random_image(shape, seed=0, dtype=np.float32): + rng = np.random.RandomState(seed) + return rng.rand(*shape).astype(dtype) + + +# --------------------------------------------------------------------------- +# Smoothing +# --------------------------------------------------------------------------- + +@pytest.mark.parametrize("shape", SHAPES_2D + SHAPES_3D) +@pytest.mark.parametrize("sigma", SIGMAS) +def test_gaussian_smoothing_matches_scipy(shape, sigma): + img = _random_image(shape) + got = bf.gaussian_smoothing(img, sigma) + ref = ndimage.gaussian_filter(img, sigma, mode="mirror") + assert got.shape == ref.shape + assert got.dtype == np.float32 + np.testing.assert_allclose(got, ref, atol=1e-3) + + +def test_gaussian_smoothing_anisotropic_sigma(): + img = _random_image((32, 32)) + got = bf.gaussian_smoothing(img, [1.0, 2.5]) + ref = ndimage.gaussian_filter(img, [1.0, 2.5], mode="mirror") + np.testing.assert_allclose(got, ref, atol=1e-3) + + +def test_gaussian_smoothing_3d_anisotropic_sigma(): + vol = _random_image((8, 12, 16)) + got = bf.gaussian_smoothing(vol, [0.7, 1.2, 2.1]) + ref = ndimage.gaussian_filter(vol, [0.7, 1.2, 2.1], mode="mirror") + np.testing.assert_allclose(got, ref, atol=1e-3) + + +# --------------------------------------------------------------------------- +# Derivative +# --------------------------------------------------------------------------- + +@pytest.mark.parametrize("order", [[1, 0], [0, 1], [2, 0], [0, 2], [1, 1]]) +def test_gaussian_derivative_2d_matches_scipy(order): + img = _random_image((32, 48)) + got = bf.gaussian_derivative(img, 1.5, order) + ref = ndimage.gaussian_filter(img, 1.5, order=order, mode="mirror") + np.testing.assert_allclose(got, ref, atol=1e-3) + + +@pytest.mark.parametrize("order", [[1, 0, 0], [0, 1, 0], [0, 0, 1], [2, 0, 0], + [1, 1, 0], [1, 0, 1]]) +def test_gaussian_derivative_3d_matches_scipy(order): + vol = _random_image((8, 12, 16)) + got = bf.gaussian_derivative(vol, 1.2, order) + ref = ndimage.gaussian_filter(vol, 1.2, order=order, mode="mirror") + np.testing.assert_allclose(got, ref, atol=1e-3) + + +# --------------------------------------------------------------------------- +# Gradient magnitude +# --------------------------------------------------------------------------- + +@pytest.mark.parametrize("shape", SHAPES_2D + SHAPES_3D) +def test_gradient_magnitude_matches_scipy(shape): + img = _random_image(shape) + got = bf.gaussian_gradient_magnitude(img, 1.5) + ref = ndimage.gaussian_gradient_magnitude(img, 1.5, mode="mirror") + np.testing.assert_allclose(got, ref, atol=1e-3) + + +# --------------------------------------------------------------------------- +# Laplacian of Gaussian +# --------------------------------------------------------------------------- + +@pytest.mark.parametrize("shape", SHAPES_2D + SHAPES_3D) +def test_laplacian_of_gaussian_matches_scipy(shape): + img = _random_image(shape) + got = bf.laplacian_of_gaussian(img, 1.5) + ref = ndimage.gaussian_laplace(img, 1.5, mode="mirror") + np.testing.assert_allclose(got, ref, atol=1e-3) + + +# --------------------------------------------------------------------------- +# Hessian eigenvalues +# --------------------------------------------------------------------------- + +def _hessian_eigenvalues_reference_2d(img, sigma): + hyy = ndimage.gaussian_filter(img, sigma, order=[2, 0], mode="mirror") + hyx = ndimage.gaussian_filter(img, sigma, order=[1, 1], mode="mirror") + hxx = ndimage.gaussian_filter(img, sigma, order=[0, 2], mode="mirror") + mat = np.stack( + [np.stack([hyy, hyx], axis=-1), np.stack([hyx, hxx], axis=-1)], + axis=-2, + ) + return np.linalg.eigvalsh(mat)[..., ::-1].astype(np.float32) + + +def _hessian_eigenvalues_reference_3d(vol, sigma): + hzz = ndimage.gaussian_filter(vol, sigma, order=[2, 0, 0], mode="mirror") + hzy = ndimage.gaussian_filter(vol, sigma, order=[1, 1, 0], mode="mirror") + hzx = ndimage.gaussian_filter(vol, sigma, order=[1, 0, 1], mode="mirror") + hyy = ndimage.gaussian_filter(vol, sigma, order=[0, 2, 0], mode="mirror") + hyx = ndimage.gaussian_filter(vol, sigma, order=[0, 1, 1], mode="mirror") + hxx = ndimage.gaussian_filter(vol, sigma, order=[0, 0, 2], mode="mirror") + mat = np.stack([ + np.stack([hzz, hzy, hzx], axis=-1), + np.stack([hzy, hyy, hyx], axis=-1), + np.stack([hzx, hyx, hxx], axis=-1), + ], axis=-2) + return np.linalg.eigvalsh(mat)[..., ::-1].astype(np.float32) + + +@pytest.mark.parametrize("shape", SHAPES_2D) +def test_hessian_eigenvalues_2d_matches_reference(shape): + img = _random_image(shape) + got = bf.hessian_of_gaussian_eigenvalues(img, 1.5) + ref = _hessian_eigenvalues_reference_2d(img, 1.5) + assert got.shape == img.shape + (2,) + np.testing.assert_allclose(got, ref, atol=2e-3) + + +@pytest.mark.parametrize("shape", SHAPES_3D) +def test_hessian_eigenvalues_3d_matches_reference(shape): + vol = _random_image(shape) + got = bf.hessian_of_gaussian_eigenvalues(vol, 1.2) + ref = _hessian_eigenvalues_reference_3d(vol, 1.2) + assert got.shape == vol.shape + (3,) + # 3x3 trig closed-form needs a slightly looser tolerance. + np.testing.assert_allclose(got, ref, atol=5e-3) + + +def test_hessian_eigenvalues_descending_order_2d(): + img = _random_image((32, 32)) + got = bf.hessian_of_gaussian_eigenvalues(img, 1.5) + assert np.all(got[..., 0] >= got[..., 1]) + + +def test_hessian_eigenvalues_descending_order_3d(): + vol = _random_image((8, 16, 16)) + got = bf.hessian_of_gaussian_eigenvalues(vol, 1.2) + assert np.all(got[..., 0] >= got[..., 1]) + assert np.all(got[..., 1] >= got[..., 2]) + + +# --------------------------------------------------------------------------- +# Structure tensor eigenvalues +# --------------------------------------------------------------------------- + +def _structure_tensor_eigenvalues_reference_2d(img, inner, outer): + gy = ndimage.gaussian_filter(img, inner, order=[1, 0], mode="mirror") + gx = ndimage.gaussian_filter(img, inner, order=[0, 1], mode="mirror") + syy = ndimage.gaussian_filter(gy * gy, outer, mode="mirror") + syx = ndimage.gaussian_filter(gy * gx, outer, mode="mirror") + sxx = ndimage.gaussian_filter(gx * gx, outer, mode="mirror") + mat = np.stack( + [np.stack([syy, syx], axis=-1), np.stack([syx, sxx], axis=-1)], + axis=-2, + ) + return np.linalg.eigvalsh(mat)[..., ::-1].astype(np.float32) + + +def test_structure_tensor_eigenvalues_2d_matches_reference(): + img = _random_image((48, 48)) + got = bf.structure_tensor_eigenvalues(img, 1.0, 2.0) + ref = _structure_tensor_eigenvalues_reference_2d(img, 1.0, 2.0) + assert got.shape == img.shape + (2,) + np.testing.assert_allclose(got, ref, atol=2e-3) + + +def test_structure_tensor_eigenvalues_3d_shape_and_order(): + vol = _random_image((8, 16, 16)) + got = bf.structure_tensor_eigenvalues(vol, 1.0, 2.0) + assert got.shape == vol.shape + (3,) + assert np.all(got[..., 0] >= got[..., 1]) + assert np.all(got[..., 1] >= got[..., 2]) + # All eigenvalues of a positive-semidefinite tensor are >= 0. + assert np.all(got >= -1e-6) + + +# --------------------------------------------------------------------------- +# dtype handling +# --------------------------------------------------------------------------- + +def test_float64_input_returns_float64(): + img = _random_image((32, 32), dtype=np.float64) + got = bf.gaussian_smoothing(img, 1.0) + assert got.dtype == np.float64 + ref = ndimage.gaussian_filter(img.astype(np.float32), 1.0, mode="mirror") + np.testing.assert_allclose(got, ref.astype(np.float64), atol=1e-3) + + +def test_uint8_input_returns_float32(): + rng = np.random.RandomState(0) + img = rng.randint(0, 256, size=(32, 32), dtype=np.uint8) + got = bf.gaussian_smoothing(img, 1.0) + assert got.dtype == np.float32 + ref = ndimage.gaussian_filter(img.astype(np.float32), 1.0, mode="mirror") + # uint8 values can be up to 255; float32 accumulation produces O(0.1) + # differences at that magnitude. We only care that the result is in the + # right ballpark. + np.testing.assert_allclose(got, ref, atol=0.1) + + +def test_uint16_input_returns_float32(): + rng = np.random.RandomState(0) + img = rng.randint(0, 4096, size=(32, 32), dtype=np.uint16) + got = bf.gaussian_smoothing(img, 1.0) + assert got.dtype == np.float32 + + +def test_non_contiguous_input_is_handled(): + img = _random_image((32, 32)) + sliced = img[::2, ::2] # non-contiguous strided view + got = bf.gaussian_smoothing(sliced, 1.0) + ref = ndimage.gaussian_filter(sliced, 1.0, mode="mirror") + np.testing.assert_allclose(got, ref, atol=1e-3) + + +# --------------------------------------------------------------------------- +# Error paths +# --------------------------------------------------------------------------- + +def test_wrong_ndim_raises(): + with pytest.raises(ValueError, match="must be 2D or 3D"): + bf.gaussian_smoothing(np.zeros(8, dtype=np.float32), 1.0) + with pytest.raises(ValueError, match="must be 2D or 3D"): + bf.gaussian_smoothing(np.zeros((4, 4, 4, 4), dtype=np.float32), 1.0) + + +def test_unsupported_dtype_raises(): + with pytest.raises(TypeError, match="dtype"): + bf.gaussian_smoothing(np.zeros((8, 8), dtype=np.int32), 1.0) + + +def test_non_positive_sigma_raises(): + img = _random_image((16, 16)) + with pytest.raises(Exception): # noqa: B017 - C++ -> invalid_argument + bf.gaussian_smoothing(img, 0.0) + with pytest.raises(Exception): # noqa: B017 + bf.gaussian_smoothing(img, -1.0) + + +def test_sigma_length_mismatch_raises(): + img = _random_image((16, 16)) + with pytest.raises(ValueError, match="sigma"): + bf.gaussian_smoothing(img, [1.0, 2.0, 3.0]) + + +def test_invalid_order_raises(): + img = _random_image((16, 16)) + with pytest.raises(Exception): # noqa: B017 - C++ -> invalid_argument + bf.gaussian_derivative(img, 1.0, 3) + with pytest.raises(Exception): # noqa: B017 + bf.gaussian_derivative(img, 1.0, -1) From f462a2a446d3b1711fd102d5b268258d69f6f135 Mon Sep 17 00:00:00 2001 From: Constantin Pape Date: Sun, 17 May 2026 17:30:02 +0200 Subject: [PATCH 2/3] Add filter benchmarks --- development/filters/PERFORMANCE_NOTES.md | 165 +++++++ development/filters/_bench_utils.py | 562 +++++++++++++++++++++++ development/filters/benchmark.py | 165 +++++++ development/filters/check_parity.py | 134 ++++++ 4 files changed, 1026 insertions(+) create mode 100644 development/filters/PERFORMANCE_NOTES.md create mode 100644 development/filters/_bench_utils.py create mode 100644 development/filters/benchmark.py create mode 100644 development/filters/check_parity.py diff --git a/development/filters/PERFORMANCE_NOTES.md b/development/filters/PERFORMANCE_NOTES.md new file mode 100644 index 0000000..063ace2 --- /dev/null +++ b/development/filters/PERFORMANCE_NOTES.md @@ -0,0 +1,165 @@ +# Filter benchmark — performance notes + +Results from `development/filters/benchmark.py` for the six Tier-1 filters +on 2D and 3D test data from `skimage.data`. Re-run with: + +```bash +python development/filters/check_parity.py # parity gate, must PASS first +python development/filters/benchmark.py # headline numbers +``` + +All four libraries produce the same response (verified by `check_parity.py` +on the image interior, dropping a `window_size * sigma` margin on each side; +or `inner + outer` for the structure tensor). Reported numbers are +**median wall-clock per call** across 5 timed repeats, after one untimed +warmup, with calls interleaved round-robin between libraries to share cache +state fairly. + +## Setup + +- CPU: Intel Core i7-1185G7 (Tiger Lake, 4C/8T, AVX2 + AVX-512) +- Compiler: gcc 15.2.0 (conda-forge), `-O3`, no `-march=native` +- Python 3.11.15 on Linux x86_64 +- `numpy 2.4.5`, `scipy 1.17.1`, `vigra 1.12.3`, `fastfilters 0.3-4-g87f08b5`, + `bioimage_cpp 0.1.0` +- All libraries called single-threaded; `bioimage_cpp` has no SIMD intrinsics + (Tier-1: scalar C++20 + compiler auto-vectorisation). + +## Benchmark configuration + +| parameter | value | notes | +|---|---|---| +| sigma | 1.5 | smoothing / derivative / LoG / gradient / Hessian | +| inner_sigma | 1.0 | structure-tensor gradient scale | +| outer_sigma | 2.0 | structure-tensor smoothing scale | +| window_size | 3.0 | kernel half-width / sigma; matched across libraries | +| repeats | 5 | timed, interleaved round-robin | +| 2D image | `skimage.data.camera()` | 512×512 → float32, normalised to [0, 1] | +| 3D volume | `skimage.data.cells3d()[:, 1]` | 60×256×256 nuclei channel → float32 | + +`gaussian_derivative` uses order = 1 along the trailing axis only. `fastfilters` +only supports uniform per-axis order, so its row is `n/a`. + +## Headline ratios + +Geometric mean of `bioimage_cpp.median / other.median` across all benched +(filter, dim) combinations. **>1.0 means `bioimage_cpp` is slower** than the +other library; **<1.0 means faster**. + +| comparison | geomean ratio | n | +|---|---|---| +| `bioimage_cpp` / `fastfilters` | **2.00** | 10 | +| `bioimage_cpp` / `vigra` | **0.18** | 12 | +| `bioimage_cpp` / `scipy` | **0.20** | 12 | + +Interpretation: `fastfilters` (hand-AVX2) is 2× faster than ours on average; +ours is ~5× faster than vigra and ~5× faster than scipy/numpy. + +## 2D results — `camera()` 512×512 + +Times in milliseconds (median of 5 repeats). The `x ours` column is +`bioimage_cpp.median / this_lib.median`; **values >1.0 mean the library is +faster than `bioimage_cpp`**. + +| filter | bioimage_cpp ms | fastfilters ms | x ours | vigra ms | x ours | scipy ms | x ours | +|---|---:|---:|---:|---:|---:|---:|---:| +| gaussian_smoothing | 0.82 | 0.49 | 1.67 | 5.18 | 0.16 | 3.59 | 0.23 | +| gaussian_derivative | 0.83 | n/a | n/a | 5.15 | 0.16 | 3.85 | 0.21 | +| gaussian_gradient_magnitude | 2.11 | 0.91 | 2.32 | 13.85 | 0.15 | 7.77 | 0.27 | +| laplacian_of_gaussian | 1.92 | 0.96 | 2.00 | 8.88 | 0.22 | 7.50 | 0.26 | +| hessian_of_gaussian_eigenvalues | 3.02 | 1.67 | 1.81 | 19.55 | 0.15 | 77.79 | 0.04 | +| structure_tensor_eigenvalues | 5.16 | 2.86 | 1.80 | 26.22 | 0.20 | 82.99 | 0.06 | + +## 3D results — `cells3d()[:, 1]` 60×256×256 + +| filter | bioimage_cpp ms | fastfilters ms | x ours | vigra ms | x ours | scipy ms | x ours | +|---|---:|---:|---:|---:|---:|---:|---:| +| gaussian_smoothing | 33.27 | 15.02 | 2.22 | 192.29 | 0.17 | 85.37 | 0.39 | +| gaussian_derivative | 33.05 | n/a | n/a | 195.27 | 0.17 | 86.19 | 0.38 | +| gaussian_gradient_magnitude | 106.42 | 61.91 | 1.72 | 640.83 | 0.17 | 276.64 | 0.38 | +| laplacian_of_gaussian | 103.69 | 61.93 | 1.67 | 584.71 | 0.18 | 266.87 | 0.39 | +| hessian_of_gaussian_eigenvalues | 486.38 | 175.75 | 2.77 | 2195.57 | 0.22 | 4031.79 | 0.12 | +| structure_tensor_eigenvalues | 599.90 | 263.14 | 2.28 | 1938.37 | 0.31 | 4301.12 | 0.14 | + +## Per-filter throughput vs `bioimage_cpp` (megapixels / second) + +Throughput = `image_size / median_time`. Useful for cross-shape comparison. + +### 2D (512×512 = 0.262 megapixels per output) + +| filter | bioimage_cpp | fastfilters | vigra | scipy | +|---|---:|---:|---:|---:| +| gaussian_smoothing | 320 | 535 | 51 | 73 | +| gaussian_gradient_magnitude | 124 | 288 | 19 | 34 | +| laplacian_of_gaussian | 136 | 273 | 30 | 35 | +| hessian_of_gaussian_eigenvalues | 87 | 157 | 13 | 3 | +| structure_tensor_eigenvalues | 51 | 92 | 10 | 3 | + +### 3D (60×256×256 = 3.93 megavoxels per output) + +| filter | bioimage_cpp | fastfilters | vigra | scipy | +|---|---:|---:|---:|---:| +| gaussian_smoothing | 118 | 262 | 20 | 46 | +| gaussian_gradient_magnitude | 37 | 64 | 6 | 14 | +| laplacian_of_gaussian | 38 | 64 | 7 | 15 | +| hessian_of_gaussian_eigenvalues | 8 | 22 | 2 | 1 | +| structure_tensor_eigenvalues | 7 | 15 | 2 | 1 | + +## Takeaways + +- **vs `fastfilters`** (the hand-AVX2 target). We sit at 1.67× – 2.77× + slower across the board, exactly in the 1.0×–2.0× band the Tier-1 plan + predicted (a touch worse on Hessian/structure-tensor 3D where eigenvalue + arithmetic dominates and benefits less from auto-vectorisation). This is + the gap a Tier-2 manual-AVX2 path would aim to close. +- **vs `vigra`**. ~5× faster across the board. Same algorithmic family, + scalar code in both cases; the win comes from fewer abstraction layers + and tighter kernels (constant-bound inner loops, half-kernel storage, + X-strip pattern for the strided pass). +- **vs `scipy.ndimage` (+ `numpy.linalg.eigvalsh` for the eigenvalue + filters)**. ~5× faster on simple filters and ~10–25× faster on the + eigenvalue paths. `scipy.ndimage` itself is competitive on plain + convolution; the eigenvalue gap is almost entirely the per-pixel numpy + `eigvalsh` cost — exactly what scipy-only users currently pay. + +## Where the Tier-2 SIMD work would help most + +Largest absolute gaps to `fastfilters` (3D, in ms per call): + +| filter | ours | ff | absolute gap | gap × 5 calls/sec | +|---|---:|---:|---:|---:| +| hessian_of_gaussian_eigenvalues | 486 | 176 | 311 ms | 1.55 s/sec saved | +| structure_tensor_eigenvalues | 600 | 263 | 337 ms | 1.69 s/sec saved | +| gaussian_gradient_magnitude | 106 | 62 | 44 ms | 0.22 s/sec saved | + +If/when Tier-2 lands, instrument the Hessian and structure-tensor 3D paths +first — they share the same separable-FIR primitives plus the 3×3 trig +eigensolver, so closing the gap on those carries the gradient/magnitude / +LoG paths along for free. + +## Reproducibility notes + +- Run after `pip install -e . --no-build-isolation` so the C++ extension + matches the source tree. +- For absolute-time comparisons across machines, also report CPU model, + compiler version, and whether `-march=native` was set (we do NOT set it + in normal builds; this benchmark used the default `-O3`). +- For a quick re-check use `--small` (crops to 128×128 and 32×64×64); + finishes in a few seconds. +- For raw per-(filter, library) timings (useful for plotting), pass + `--csv path.csv`. + +## Known caveats reflected in the adapters + +- `fastfilters.gaussianDerivative` only accepts a uniform per-axis order; + the bench cell is `n/a` rather than running an unequal operation. +- `fastfilters.structureTensorEigenvalues` swaps `innerScale` / + `outerScale` at the Python boundary (its wrapper calls the C function + with the args in the opposite order — `src/python/core.cxx:328` vs + `src/library/fir_filters.c:156` in the fastfilters source). The adapter + swaps them back so the bench compares the same operation as `vigra` / + `bioimage_cpp`. +- `scipy` and `bioimage_cpp` use scipy-style `mirror` (reflect without + edge-pixel repeat); `vigra` / `fastfilters` use reflect with edge + repeat. Parity is checked on the image interior to absorb the + difference. diff --git a/development/filters/_bench_utils.py b/development/filters/_bench_utils.py new file mode 100644 index 0000000..24c1d12 --- /dev/null +++ b/development/filters/_bench_utils.py @@ -0,0 +1,562 @@ +"""Shared helpers for the filters benchmark scripts. + +This module is intentionally not part of the test suite. It provides: + +* data loaders that prepare the same float32 array for every library +* one adapter per (library, filter) that hides parameter-name and output-shape + differences so the harness can call any adapter as ``fn(image)`` +* a small interleaved timing harness and a fixed-width report formatter +""" + +from __future__ import annotations + +from collections.abc import Callable, Sequence +from dataclasses import dataclass +from statistics import median +from time import perf_counter + +import numpy as np + + +LIBRARIES: tuple[str, ...] = ("bioimage_cpp", "fastfilters", "vigra", "scipy") +FILTERS: tuple[str, ...] = ( + "gaussian_smoothing", + "gaussian_derivative", + "gaussian_gradient_magnitude", + "laplacian_of_gaussian", + "hessian_of_gaussian_eigenvalues", + "structure_tensor_eigenvalues", +) + +# Per-axis order tuple used for gaussian_derivative. Convention: first +# derivative along the trailing spatial axis ("d/dx"), zero elsewhere. +DERIVATIVE_ORDER_2D = (0, 1) +DERIVATIVE_ORDER_3D = (0, 0, 1) + +NOT_APPLICABLE = "n/a" + + +# --------------------------------------------------------------------------- +# Test data +# --------------------------------------------------------------------------- + +def load_2d(*, crop: tuple[int, int] | None = None) -> np.ndarray: + """skimage.data.camera() as contiguous float32 in [0, 1].""" + from skimage import data + arr = data.camera() + if crop is not None: + arr = arr[: crop[0], : crop[1]] + return np.ascontiguousarray(arr.astype(np.float32) / 255.0) + + +def load_3d(*, crop: tuple[int, int, int] | None = None) -> np.ndarray: + """skimage.data.cells3d() nuclei channel as contiguous float32 in [0, 1].""" + from skimage import data + vol = data.cells3d()[:, 1] # (60, 256, 256), uint16 + if crop is not None: + vol = vol[: crop[0], : crop[1], : crop[2]] + arr = vol.astype(np.float32) + arr /= float(arr.max() if arr.max() > 0 else 1.0) + return np.ascontiguousarray(arr) + + +# --------------------------------------------------------------------------- +# Bench config: knobs every adapter must respect for fair comparison. +# --------------------------------------------------------------------------- + +@dataclass(frozen=True) +class BenchConfig: + sigma: float = 1.5 + inner_sigma: float = 1.0 + outer_sigma: float = 2.0 + # window_size in vigra / fastfilters / bioimage_cpp and truncate in scipy + # both mean "kernel half-width / sigma". Fixed to make kernel sizes match + # across libraries. + window_size: float = 3.0 + + @property + def truncate(self) -> float: + return self.window_size + + def derivative_order(self, ndim: int) -> tuple[int, ...]: + return DERIVATIVE_ORDER_2D if ndim == 2 else DERIVATIVE_ORDER_3D + + +# --------------------------------------------------------------------------- +# Adapters: each builder returns a callable fn(image) -> np.ndarray. +# A None return value indicates the library does not support that filter. +# --------------------------------------------------------------------------- + +# ---- bioimage_cpp adapters ------------------------------------------------ + +def _bic_smoothing(cfg: BenchConfig): + from bioimage_cpp import filters as bf + sigma, ws = cfg.sigma, cfg.window_size + def fn(image): + return bf.gaussian_smoothing(image, sigma, window_size=ws) + return fn + + +def _bic_derivative(cfg: BenchConfig): + from bioimage_cpp import filters as bf + sigma, ws = cfg.sigma, cfg.window_size + def fn(image): + order = cfg.derivative_order(image.ndim) + return bf.gaussian_derivative(image, sigma, order, window_size=ws) + return fn + + +def _bic_grad_magnitude(cfg: BenchConfig): + from bioimage_cpp import filters as bf + sigma, ws = cfg.sigma, cfg.window_size + def fn(image): + return bf.gaussian_gradient_magnitude(image, sigma, window_size=ws) + return fn + + +def _bic_log(cfg: BenchConfig): + from bioimage_cpp import filters as bf + sigma, ws = cfg.sigma, cfg.window_size + def fn(image): + return bf.laplacian_of_gaussian(image, sigma, window_size=ws) + return fn + + +def _bic_hessian_ev(cfg: BenchConfig): + from bioimage_cpp import filters as bf + sigma, ws = cfg.sigma, cfg.window_size + def fn(image): + return bf.hessian_of_gaussian_eigenvalues(image, sigma, window_size=ws) + return fn + + +def _bic_st_ev(cfg: BenchConfig): + from bioimage_cpp import filters as bf + inner, outer, ws = cfg.inner_sigma, cfg.outer_sigma, cfg.window_size + def fn(image): + return bf.structure_tensor_eigenvalues(image, inner, outer, window_size=ws) + return fn + + +# ---- fastfilters adapters ------------------------------------------------- + +def _ff_smoothing(cfg: BenchConfig): + import fastfilters as ff + sigma, ws = cfg.sigma, cfg.window_size + def fn(image): + return ff.gaussianSmoothing(image, sigma, window_size=ws) + return fn + + +def _ff_derivative(cfg: BenchConfig): + # fastfilters only supports uniform per-axis order; our derivative + # benchmark is order=1 along the last axis only, so fastfilters cannot + # match the operation. Mark as not applicable. + return None + + +def _ff_grad_magnitude(cfg: BenchConfig): + import fastfilters as ff + sigma, ws = cfg.sigma, cfg.window_size + def fn(image): + return ff.gaussianGradientMagnitude(image, sigma, window_size=ws) + return fn + + +def _ff_log(cfg: BenchConfig): + import fastfilters as ff + sigma, ws = cfg.sigma, cfg.window_size + def fn(image): + return ff.laplacianOfGaussian(image, sigma, window_size=ws) + return fn + + +def _ff_hessian_ev(cfg: BenchConfig): + import fastfilters as ff + sigma, ws = cfg.sigma, cfg.window_size + def fn(image): + return ff.hessianOfGaussianEigenvalues(image, sigma, window_size=ws) + return fn + + +def _ff_st_ev(cfg: BenchConfig): + import fastfilters as ff + inner, outer, ws = cfg.inner_sigma, cfg.outer_sigma, cfg.window_size + # fastfilters' Python wrapper calls the C function as + # fastfilters_fir_structure_tensor2d(in, sigma_inner, sigma_outer, ...) + # but the C signature is (sigma_outer, sigma_inner) — see src/python/core.cxx + # line 328 vs src/library/fir_filters.c line 156 in the fastfilters source. + # The two scales are therefore swapped at the Python boundary relative to + # vigra / scipy / bioimage_cpp. Swap them back so this adapter computes + # the same operation as the rest. + def fn(image): + return ff.structureTensorEigenvalues(image, innerScale=outer, outerScale=inner, window_size=ws) + return fn + + +# ---- vigra adapters ------------------------------------------------------- + +def _vigra_smoothing(cfg: BenchConfig): + import vigra.filters as vf + sigma, ws = cfg.sigma, cfg.window_size + def fn(image): + return np.asarray(vf.gaussianSmoothing(image, sigma, window_size=ws)) + return fn + + +def _vigra_derivative(cfg: BenchConfig): + import vigra.filters as vf + sigma, ws = cfg.sigma, cfg.window_size + def fn(image): + order = list(cfg.derivative_order(image.ndim)) + return np.asarray(vf.gaussianDerivative(image, sigma, order, window_size=ws)) + return fn + + +def _vigra_grad_magnitude(cfg: BenchConfig): + import vigra.filters as vf + sigma, ws = cfg.sigma, cfg.window_size + def fn(image): + return np.asarray(vf.gaussianGradientMagnitude(image, sigma, window_size=ws)) + return fn + + +def _vigra_log(cfg: BenchConfig): + import vigra.filters as vf + sigma, ws = cfg.sigma, cfg.window_size + def fn(image): + return np.asarray(vf.laplacianOfGaussian(image, scale=sigma, window_size=ws)) + return fn + + +def _vigra_hessian_ev(cfg: BenchConfig): + import vigra.filters as vf + sigma, ws = cfg.sigma, cfg.window_size + def fn(image): + return np.asarray(vf.hessianOfGaussianEigenvalues(image, scale=sigma, window_size=ws)) + return fn + + +def _vigra_st_ev(cfg: BenchConfig): + import vigra.filters as vf + inner, outer, ws = cfg.inner_sigma, cfg.outer_sigma, cfg.window_size + def fn(image): + return np.asarray( + vf.structureTensorEigenvalues( + image, innerScale=inner, outerScale=outer, window_size=ws, + ) + ) + return fn + + +# ---- scipy / numpy adapters ---------------------------------------------- +# +# The two eigenvalue filters do not have a direct scipy entry point; we hand- +# build the components using scipy.ndimage gaussian filters with matching +# kernel support and use numpy.linalg.eigvalsh per pixel. This is what a +# scipy user would write and represents the realistic scipy-only baseline. + +def _scipy_smoothing(cfg: BenchConfig): + from scipy import ndimage + sigma, tr = cfg.sigma, cfg.truncate + def fn(image): + return ndimage.gaussian_filter(image, sigma, mode="mirror", truncate=tr) + return fn + + +def _scipy_derivative(cfg: BenchConfig): + from scipy import ndimage + sigma, tr = cfg.sigma, cfg.truncate + def fn(image): + order = list(cfg.derivative_order(image.ndim)) + return ndimage.gaussian_filter(image, sigma, order=order, mode="mirror", truncate=tr) + return fn + + +def _scipy_grad_magnitude(cfg: BenchConfig): + from scipy import ndimage + sigma, tr = cfg.sigma, cfg.truncate + def fn(image): + return ndimage.gaussian_gradient_magnitude( + image, sigma, mode="mirror", truncate=tr, + ) + return fn + + +def _scipy_log(cfg: BenchConfig): + from scipy import ndimage + sigma, tr = cfg.sigma, cfg.truncate + def fn(image): + return ndimage.gaussian_laplace(image, sigma, mode="mirror", truncate=tr) + return fn + + +def _scipy_hessian_ev(cfg: BenchConfig): + from scipy import ndimage + sigma, tr = cfg.sigma, cfg.truncate + + def _components_2d(image): + hyy = ndimage.gaussian_filter(image, sigma, order=[2, 0], mode="mirror", truncate=tr) + hyx = ndimage.gaussian_filter(image, sigma, order=[1, 1], mode="mirror", truncate=tr) + hxx = ndimage.gaussian_filter(image, sigma, order=[0, 2], mode="mirror", truncate=tr) + return hyy, hyx, hxx + + def _components_3d(image): + return ( + ndimage.gaussian_filter(image, sigma, order=[2, 0, 0], mode="mirror", truncate=tr), + ndimage.gaussian_filter(image, sigma, order=[1, 1, 0], mode="mirror", truncate=tr), + ndimage.gaussian_filter(image, sigma, order=[1, 0, 1], mode="mirror", truncate=tr), + ndimage.gaussian_filter(image, sigma, order=[0, 2, 0], mode="mirror", truncate=tr), + ndimage.gaussian_filter(image, sigma, order=[0, 1, 1], mode="mirror", truncate=tr), + ndimage.gaussian_filter(image, sigma, order=[0, 0, 2], mode="mirror", truncate=tr), + ) + + def fn(image): + if image.ndim == 2: + hyy, hyx, hxx = _components_2d(image) + mat = np.stack( + [np.stack([hyy, hyx], axis=-1), np.stack([hyx, hxx], axis=-1)], + axis=-2, + ) + else: + hzz, hzy, hzx, hyy, hyx, hxx = _components_3d(image) + mat = np.stack([ + np.stack([hzz, hzy, hzx], axis=-1), + np.stack([hzy, hyy, hyx], axis=-1), + np.stack([hzx, hyx, hxx], axis=-1), + ], axis=-2) + evs = np.linalg.eigvalsh(mat)[..., ::-1] + return evs.astype(np.float32, copy=False) + + return fn + + +def _scipy_st_ev(cfg: BenchConfig): + from scipy import ndimage + inner, outer, tr = cfg.inner_sigma, cfg.outer_sigma, cfg.truncate + + def _grads(image): + unit = np.eye(image.ndim, dtype=int) + return [ + ndimage.gaussian_filter(image, inner, order=list(unit[i]), + mode="mirror", truncate=tr) + for i in range(image.ndim) + ] + + def fn(image): + grads = _grads(image) + comps = {} + for i in range(image.ndim): + for j in range(i, image.ndim): + comps[(i, j)] = ndimage.gaussian_filter( + grads[i] * grads[j], outer, mode="mirror", truncate=tr, + ) + if image.ndim == 2: + mat = np.stack([ + np.stack([comps[(0, 0)], comps[(0, 1)]], axis=-1), + np.stack([comps[(0, 1)], comps[(1, 1)]], axis=-1), + ], axis=-2) + else: + mat = np.stack([ + np.stack([comps[(0, 0)], comps[(0, 1)], comps[(0, 2)]], axis=-1), + np.stack([comps[(0, 1)], comps[(1, 1)], comps[(1, 2)]], axis=-1), + np.stack([comps[(0, 2)], comps[(1, 2)], comps[(2, 2)]], axis=-1), + ], axis=-2) + evs = np.linalg.eigvalsh(mat)[..., ::-1] + return evs.astype(np.float32, copy=False) + + return fn + + +# --------------------------------------------------------------------------- +# Adapter table +# --------------------------------------------------------------------------- + +ADAPTERS: dict[str, dict[str, Callable[[BenchConfig], Callable[[np.ndarray], np.ndarray] | None]]] = { + "gaussian_smoothing": { + "bioimage_cpp": _bic_smoothing, + "fastfilters": _ff_smoothing, + "vigra": _vigra_smoothing, + "scipy": _scipy_smoothing, + }, + "gaussian_derivative": { + "bioimage_cpp": _bic_derivative, + "fastfilters": _ff_derivative, + "vigra": _vigra_derivative, + "scipy": _scipy_derivative, + }, + "gaussian_gradient_magnitude": { + "bioimage_cpp": _bic_grad_magnitude, + "fastfilters": _ff_grad_magnitude, + "vigra": _vigra_grad_magnitude, + "scipy": _scipy_grad_magnitude, + }, + "laplacian_of_gaussian": { + "bioimage_cpp": _bic_log, + "fastfilters": _ff_log, + "vigra": _vigra_log, + "scipy": _scipy_log, + }, + "hessian_of_gaussian_eigenvalues": { + "bioimage_cpp": _bic_hessian_ev, + "fastfilters": _ff_hessian_ev, + "vigra": _vigra_hessian_ev, + "scipy": _scipy_hessian_ev, + }, + "structure_tensor_eigenvalues": { + "bioimage_cpp": _bic_st_ev, + "fastfilters": _ff_st_ev, + "vigra": _vigra_st_ev, + "scipy": _scipy_st_ev, + }, +} + + +def build_adapters(filter_name: str, cfg: BenchConfig) -> dict[str, Callable | None]: + """Build a dict {library: callable_or_None} for one filter.""" + return {lib: builder(cfg) for lib, builder in ADAPTERS[filter_name].items()} + + +# --------------------------------------------------------------------------- +# Timing +# --------------------------------------------------------------------------- + +def time_interleaved( + callables: dict[str, Callable[[np.ndarray], np.ndarray]], + image: np.ndarray, + repeats: int, +) -> dict[str, dict]: + """Run each callable ``repeats`` times in round-robin order. + + Returns {library: {"timings": [s, ...], "median": s, "min": s, + "result": last_result_array}}. + The order in which callables are timed is rotated by repeat index so that + no library is systematically advantaged by cache state. + """ + libs = list(callables.keys()) + # One untimed warmup call per library covers lazy init. + for fn in callables.values(): + fn(image) + + timings: dict[str, list[float]] = {lib: [] for lib in libs} + last_result: dict[str, np.ndarray] = {} + + n = len(libs) + for r in range(repeats): + rotation = r % n + order = libs[rotation:] + libs[:rotation] + for lib in order: + fn = callables[lib] + t0 = perf_counter() + result = fn(image) + t1 = perf_counter() + timings[lib].append(t1 - t0) + last_result[lib] = np.asarray(result) + + return { + lib: { + "timings": timings[lib], + "median": median(timings[lib]), + "min": min(timings[lib]), + "result": last_result[lib], + } + for lib in libs + } + + +# --------------------------------------------------------------------------- +# Reporting +# --------------------------------------------------------------------------- + +def format_results_table( + rows: list[dict], + *, + reference_library: str = "bioimage_cpp", + libraries: Sequence[str] = LIBRARIES, +) -> str: + """Render a fixed-width text table. + + Each row is a dict with at least: + - "filter": str + - "dim": str (e.g. "2D" or "3D") + - "shape": str (e.g. "(512, 512)") + - "results": dict {lib: {"median": s, "min": s} or None} + """ + headers = ["filter", "dim", "shape"] + for lib in libraries: + headers.append(f"{lib} ms") + # "x ours" reads as "speedup factor over bioimage_cpp": ours_time / + # this_lib_time, so >1 means this lib is faster than ours. + headers.append("x ours") + + str_rows: list[list[str]] = [] + for row in rows: + ref = row["results"].get(reference_library) or {} + ref_median = ref.get("median") if ref else None + line = [row["filter"], row["dim"], row["shape"]] + for lib in libraries: + r = row["results"].get(lib) + if r is None: + line.append(NOT_APPLICABLE) + line.append(NOT_APPLICABLE) + else: + line.append(f"{r['median'] * 1e3:.2f}") + if ref_median is None or lib == reference_library: + line.append("1.00" if lib == reference_library else "-") + else: + speedup = ref_median / r["median"] + line.append(f"{speedup:.2f}") + str_rows.append(line) + + # Column widths + widths = [len(h) for h in headers] + for r in str_rows: + for i, cell in enumerate(r): + widths[i] = max(widths[i], len(cell)) + + def render_row(values): + return " ".join(v.ljust(widths[i]) for i, v in enumerate(values)) + + out = [render_row(headers), render_row(["-" * w for w in widths])] + out.extend(render_row(r) for r in str_rows) + return "\n".join(out) + + +def interior_slice(image_shape: tuple[int, ...], border: int) -> tuple[slice, ...]: + """Slice that drops `border` pixels on each side of every spatial axis. + If the array has more axes than `len(image_shape)` (e.g. trailing + eigenvalue axis), those are passed through unsliced.""" + base = tuple(slice(border, size - border) for size in image_shape) + return base + + +def parity_atol_for_filter(filter_name: str) -> float: + # Per-filter float32 tolerance budgets: + # - smoothing / derivative / gradient_magnitude: dominated by single + # separable-Gaussian rounding; 2e-3 is comfortable. + # - LoG: sum of two second-derivatives accumulates a touch more noise; + # scipy's slightly different truncation rule pushes us over 2e-3. + # - eigenvalues: sqrt/acos/trig + sort tie-breaking near plateaus, + # plus two convolutions for the structure tensor. + if filter_name == "laplacian_of_gaussian": + return 5e-3 + if filter_name.endswith("_eigenvalues"): + return 5e-3 + return 2e-3 + + +def parity_border_for_filter(filter_name: str, cfg: "BenchConfig") -> int: + """How many pixels to drop on each side before comparing. + + The structure-tensor pipeline applies *two* Gaussians (gradient at + ``inner_sigma``, smoothing at ``outer_sigma``), so its boundary footprint + is the sum of both radii. Other filters only see ``sigma`` once. + """ + ws = cfg.window_size + if filter_name == "structure_tensor_eigenvalues": + return int(math.ceil(ws * (cfg.inner_sigma + cfg.outer_sigma))) + return int(math.ceil(ws * cfg.sigma)) + + +# Imports kept at the bottom for `parity_border_for_filter` annotation. +import math # noqa: E402 diff --git a/development/filters/benchmark.py b/development/filters/benchmark.py new file mode 100644 index 0000000..2b5416a --- /dev/null +++ b/development/filters/benchmark.py @@ -0,0 +1,165 @@ +"""Throughput benchmark of bioimage_cpp.filters against fastfilters, vigra, +and scipy.ndimage on 2D and 3D test data from ``skimage.data``. + +Run:: + + python development/filters/check_parity.py # gate: run me first + python development/filters/benchmark.py [--small] [--sigma 1.5] [--repeats 5] +""" + +from __future__ import annotations + +import argparse +import csv +import math +import sys +from statistics import geometric_mean + +from _bench_utils import ( + BenchConfig, + FILTERS, + LIBRARIES, + build_adapters, + format_results_table, + load_2d, + load_3d, + time_interleaved, +) + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser( + description="Benchmark bioimage_cpp.filters vs fastfilters / vigra / scipy." + ) + parser.add_argument("--sigma", type=float, default=1.5) + parser.add_argument("--inner-sigma", type=float, default=1.0) + parser.add_argument("--outer-sigma", type=float, default=2.0) + parser.add_argument("--window-size", type=float, default=3.0) + parser.add_argument("--repeats", type=int, default=5) + parser.add_argument("--small", action="store_true", + help="Crop to fast sizes for a smoke run.") + parser.add_argument("--no-3d", action="store_true") + parser.add_argument("--no-2d", action="store_true") + parser.add_argument( + "--filters", default=",".join(FILTERS), + help="Comma-separated subset of filters to benchmark.", + ) + parser.add_argument( + "--csv", default=None, + help="Optional path to write per-(filter, dim, library) median+min " + "timings as CSV.", + ) + return parser.parse_args() + + +def _load_test_data(args) -> list[tuple[str, "np.ndarray"]]: + crop_2d = (128, 128) if args.small else None + crop_3d = (32, 64, 64) if args.small else None + targets = [] + if not args.no_2d: + targets.append(("2D", load_2d(crop=crop_2d))) + if not args.no_3d: + targets.append(("3D", load_3d(crop=crop_3d))) + return targets + + +def main() -> int: + args = parse_args() + cfg = BenchConfig( + sigma=args.sigma, + inner_sigma=args.inner_sigma, + outer_sigma=args.outer_sigma, + window_size=args.window_size, + ) + requested = [f.strip() for f in args.filters.split(",") if f.strip()] + unknown = [f for f in requested if f not in FILTERS] + if unknown: + print(f"unknown filter(s): {unknown}", file=sys.stderr) + return 2 + + targets = _load_test_data(args) + + print( + f"sigma={cfg.sigma}, inner={cfg.inner_sigma}, outer={cfg.outer_sigma}, " + f"window_size={cfg.window_size}, repeats={args.repeats}" + ) + + rows = [] + csv_rows = [] + for dim_label, image in targets: + for filter_name in requested: + adapters = build_adapters(filter_name, cfg) + timed = {lib: fn for lib, fn in adapters.items() if fn is not None} + if not timed: + continue + results = time_interleaved(timed, image, repeats=args.repeats) + full_results = {lib: results.get(lib) for lib in LIBRARIES} + row = { + "filter": filter_name, + "dim": dim_label, + "shape": str(tuple(image.shape)), + "results": full_results, + } + rows.append(row) + for lib, r in full_results.items(): + if r is None: + continue + csv_rows.append({ + "filter": filter_name, + "dim": dim_label, + "shape": tuple(image.shape), + "library": lib, + "median_s": r["median"], + "min_s": r["min"], + "repeats": args.repeats, + }) + + print() + print(format_results_table(rows)) + + # Headline ratios across all benched (filter, dim) combinations. + print() + _print_headline_ratios(rows) + + if args.csv is not None: + with open(args.csv, "w", newline="") as fh: + writer = csv.DictWriter( + fh, fieldnames=["filter", "dim", "shape", "library", "median_s", "min_s", "repeats"] + ) + writer.writeheader() + for r in csv_rows: + writer.writerow(r) + print(f"wrote {args.csv}") + + return 0 + + +def _print_headline_ratios(rows: list[dict]) -> None: + """Print geometric-mean wall-time ratios. + + Each ratio is ``bioimage_cpp.median / other.median``: values above 1.0 + mean we are slower than the other library; values below 1.0 mean we are + faster. + """ + + def gm_ratio(other: str) -> tuple[float | None, int]: + ratios = [] + for row in rows: + ours = row["results"].get("bioimage_cpp") + them = row["results"].get(other) + if ours and them and them["median"] > 0: + ratios.append(ours["median"] / them["median"]) + return (geometric_mean(ratios) if ratios else None), len(ratios) + + others = ["fastfilters", "vigra", "scipy"] + print("speedup summary (geomean of bioimage_cpp.median / other.median; " + ">1.0 means bioimage_cpp slower, <1.0 means faster):") + for other in others: + gm, n = gm_ratio(other) + if gm is None: + continue + print(f" bioimage_cpp / {other:<12s} geomean = {gm:.3f} (n={n})") + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/development/filters/check_parity.py b/development/filters/check_parity.py new file mode 100644 index 0000000..2a9450f --- /dev/null +++ b/development/filters/check_parity.py @@ -0,0 +1,134 @@ +"""Equivalence check across scipy / vigra / fastfilters / bioimage_cpp filters. + +Compares every library's output against ``bioimage_cpp`` on the interior of +the image (a margin of ``2 * ceil(window_size * sigma)`` pixels is dropped on +every axis to avoid boundary-mode differences between libraries). + +Run:: + + python development/filters/check_parity.py [--sigma 1.5] [--no-3d] + +Exits non-zero on any tolerance failure. +""" + +from __future__ import annotations + +import argparse +import math +import sys + +import numpy as np + +from _bench_utils import ( + ADAPTERS, + BenchConfig, + FILTERS, + LIBRARIES, + build_adapters, + interior_slice, + load_2d, + load_3d, + parity_atol_for_filter, + parity_border_for_filter, +) + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser( + description="Equivalence check across filter implementations." + ) + parser.add_argument("--sigma", type=float, default=1.5) + parser.add_argument("--inner-sigma", type=float, default=1.0) + parser.add_argument("--outer-sigma", type=float, default=2.0) + parser.add_argument("--window-size", type=float, default=3.0) + parser.add_argument("--atol", type=float, default=None, + help="Override per-filter tolerance.") + parser.add_argument("--no-3d", action="store_true") + parser.add_argument("--no-2d", action="store_true") + parser.add_argument( + "--filters", default=",".join(FILTERS), + help="Comma-separated subset of filters to check.", + ) + return parser.parse_args() + + +def _max_abs_diff(a: np.ndarray, b: np.ndarray) -> float: + return float(np.max(np.abs(a.astype(np.float64) - b.astype(np.float64)))) + + +def _run_one(filter_name: str, image: np.ndarray, cfg: BenchConfig, atol: float) -> bool: + adapters = build_adapters(filter_name, cfg) + border = parity_border_for_filter(filter_name, cfg) + spatial = image.shape + spatial_slice = interior_slice(spatial, border) + + ref_fn = adapters["bioimage_cpp"] + ref = np.asarray(ref_fn(image)) + # The slice for the trailing eigenvalue axis (if any) is full. + full_slice = spatial_slice + (slice(None),) * (ref.ndim - len(spatial)) + ref_interior = ref[full_slice] + + all_ok = True + for lib in LIBRARIES: + if lib == "bioimage_cpp": + continue + fn = adapters[lib] + if fn is None: + print(f" {filter_name:<32s} {lib:<12s} SKIP (not supported)") + continue + got = np.asarray(fn(image)) + if got.shape != ref.shape: + print( + f" {filter_name:<32s} {lib:<12s} FAIL " + f"(shape mismatch: {got.shape} vs {ref.shape})" + ) + all_ok = False + continue + got_interior = got[full_slice] + diff = _max_abs_diff(ref_interior, got_interior) + ok = diff <= atol + status = "PASS" if ok else "FAIL" + print( + f" {filter_name:<32s} {lib:<12s} {status} max|diff|={diff:.3e}" + f" (atol={atol:.1e}, border={border})" + ) + if not ok: + all_ok = False + return all_ok + + +def main() -> int: + args = parse_args() + cfg = BenchConfig( + sigma=args.sigma, + inner_sigma=args.inner_sigma, + outer_sigma=args.outer_sigma, + window_size=args.window_size, + ) + requested = [f.strip() for f in args.filters.split(",") if f.strip()] + unknown = [f for f in requested if f not in FILTERS] + if unknown: + print(f"unknown filter(s): {unknown}", file=sys.stderr) + return 2 + + targets = [] + if not args.no_2d: + targets.append(("2D", load_2d())) + if not args.no_3d: + targets.append(("3D", load_3d())) + + any_failure = False + for dim_label, image in targets: + print(f"\n== {dim_label} parity (shape={image.shape}, dtype={image.dtype}) ==") + for filter_name in requested: + atol = args.atol if args.atol is not None else parity_atol_for_filter(filter_name) + ok = _run_one(filter_name, image, cfg, atol) + if not ok: + any_failure = True + + print("\n" + ("FAILURE" if any_failure else "OK")) + return 1 if any_failure else 0 + + +if __name__ == "__main__": + sys.exit(main()) From 090463b894ec97a842dc6ab222c09fb4c4dc8110 Mon Sep 17 00:00:00 2001 From: Constantin Pape Date: Sun, 17 May 2026 17:44:20 +0200 Subject: [PATCH 3/3] Finalize filter implementation --- development/filters/PERFORMANCE_NOTES.md | 347 +++++++++++++++++++++++ pyproject.toml | 2 +- 2 files changed, 348 insertions(+), 1 deletion(-) diff --git a/development/filters/PERFORMANCE_NOTES.md b/development/filters/PERFORMANCE_NOTES.md index 063ace2..6dddeb4 100644 --- a/development/filters/PERFORMANCE_NOTES.md +++ b/development/filters/PERFORMANCE_NOTES.md @@ -149,6 +149,353 @@ LoG paths along for free. - For raw per-(filter, library) timings (useful for plotting), pass `--csv path.csv`. +## Tier 2 SIMD — design notes (deferred) + +**Status (2026-05-17): not pursued for now.** Tier 1 sits within ~2× of +fastfilters' hand-AVX2 on the headline benchmark while being ~5× faster +than `vigra` and `scipy.ndimage`. The marginal user value of closing that +2× gap doesn't yet justify the extra build complexity and dual-path +maintenance burden. This section captures the design so a future coding +agent (or future-us) can pick it up without re-deriving the choices. + +### Trigger conditions — when to revisit + +Open this section again when **at least one** is true: + +1. Real users are hitting the Hessian-3D / structure-tensor-3D paths on + volumes large enough that ~300 ms vs ~150 ms per call materially + matters in their pipeline (typically batch feature extraction over + many large 3D blocks). +2. "Performance parity with fastfilters" becomes a stated project goal + (e.g. for a migration story or comparison documentation). +3. Profiling on a real downstream workflow shows `bioimage_cpp.filters` + is the bottleneck and the gap to fastfilters is the dominant slice. + +If none of those is true: stay on Tier 1. + +### Scope — what to ship, what to keep out + +**In scope** (the whole Tier 2 delivery): + +- Hand-written AVX2 + FMA implementations of exactly two inner kernels: + - `convolve_x_radius` — the X (innermost contiguous) pass. + - `convolve_strided_radius` — the Y/Z (strided) pass. + - Both currently live in + `include/bioimage_cpp/filters/convolve.hxx::detail`. +- One-time CPUID dispatch at module load that picks scalar vs AVX2 + function pointers for those two kernels. + +**Out of scope** (do NOT add any of these as part of Tier 2): + +- AVX-512 path. The win over AVX2 is small on memory-bound separable FIR + and doubles the kernel binary footprint; revisit only if a user with a + Sapphire Rapids / Zen 5 workload asks specifically. +- NEON / arm64 hand-tuning. Tier 1 auto-vectorization on Apple Clang is + already competitive; this would be a separate project with its own + trigger conditions. +- Replacing `std::acos` / `std::cos` in `eigenvalues.hxx` with a + vectorized math library (this is what `fastfilters` vendors as + `avx_mathfun.h`, 924 lines). The Tier-1 plan explicitly rejected + vendoring it; revisit only if eigenvalue profiling shows the trig + calls dominate the remaining gap. Don't bundle this into Tier 2. +- Any change to `kernel.hxx`, `eigenvalues.hxx`, `gaussian.hxx`, the + binding layer, or the Python wrapper. Tier 2 is a *drop-in* speedup of + two leaf functions; if you find yourself changing anything else, + something is wrong. + +### File layout + +``` +include/bioimage_cpp/filters/ + convolve.hxx # existing scalar; renamed entry points + # to point at function pointers (see below) + convolve_dispatch.hxx # NEW — function-pointer table + CPUID + +src/cpp/filters/ + convolve_avx2.cxx # NEW — AVX2+FMA kernels (compiled with + # per-file -mavx2 -mfma / /arch:AVX2) + convolve_dispatch.cxx # NEW — one-time init of the pointers +``` + +The existing `convolve_axis_x` / `convolve_axis_strided` entry points in +`convolve.hxx` keep their signatures. Their bodies switch from "directly +call `detail::convolve_x_radius`" to "call +`bioimage_cpp::filters::dispatch::convolve_x_table[R][Sym]`". Higher +levels (`gaussian.hxx`, the six composite filters, the binding layer) are +unchanged. + +### CMake wiring + +Add to the `nanobind_add_module(_core ...)` source list: + +```cmake +src/cpp/filters/convolve_avx2.cxx +src/cpp/filters/convolve_dispatch.cxx +``` + +Then attach per-file flags so only the AVX2 TU gets AVX2 instructions +(the rest of the wheel stays at the manylinux SSE2 baseline): + +```cmake +if(MSVC) + set_source_files_properties( + src/cpp/filters/convolve_avx2.cxx + PROPERTIES COMPILE_OPTIONS "/arch:AVX2" + ) +else() + set_source_files_properties( + src/cpp/filters/convolve_avx2.cxx + PROPERTIES COMPILE_OPTIONS "-mavx2;-mfma" + ) +endif() +``` + +Do **not** add `-march=native` or change the global `-O3`. The wheel +must keep installing on any pre-Haswell x86_64 machine that +manylinux2014 supports; the AVX2 instructions only execute behind the +CPUID check. + +### Runtime dispatch pattern + +In `convolve_dispatch.hxx`: + +```cpp +namespace bioimage_cpp::filters::dispatch { + +using ConvolveXFn = void (*)( + const float*, float*, std::ptrdiff_t, std::ptrdiff_t, const float* +); +using ConvolveStridedFn = void (*)( + const float*, float*, std::ptrdiff_t, std::ptrdiff_t, std::ptrdiff_t, + const float* +); + +// One entry per (radius R in 1..kMaxSpecialisedRadius, Symmetric in {0,1}). +// Filled at module load by init(). +extern ConvolveXFn convolve_x_table[kMaxSpecialisedRadius + 1][2]; +extern ConvolveStridedFn convolve_strided_table[kMaxSpecialisedRadius + 1][2]; + +void init(); // called once from bind_filters() + +} +``` + +In `convolve_dispatch.cxx`: + +```cpp +namespace { +bool detect_avx2_fma() { +#if defined(__GNUC__) || defined(__clang__) + __builtin_cpu_init(); + return __builtin_cpu_supports("avx2") && __builtin_cpu_supports("fma"); +#elif defined(_MSC_VER) + int regs1[4]; __cpuid(regs1, 1); + const bool fma = (regs1[2] & (1 << 12)) != 0; + int regs7[4]; __cpuidex(regs7, 7, 0); + const bool avx2 = (regs7[1] & (1 << 5)) != 0; + // Also OSXSAVE + XGETBV to confirm OS-saved YMM state. + ... + return avx2 && fma; +#else + return false; +#endif +} +} + +void init() { + const bool use_avx2 = detect_avx2_fma(); + // Macros generate the per-R table entries to avoid 24 hand-written + // lines (the same boost-preprocessor-style explosion fastfilters + // does, kept tiny with simple X-macros). + #define BIO_FILL(R) \ + if (use_avx2) { \ + convolve_x_table[R][0] = &avx2::convolve_x_radius_sym; \ + convolve_x_table[R][1] = &avx2::convolve_x_radius_anti; \ + convolve_strided_table[R][0] = &avx2::convolve_strided_radius_sym; \ + convolve_strided_table[R][1] = &avx2::convolve_strided_radius_anti; \ + } else { \ + convolve_x_table[R][0] = &detail::convolve_x_radius_sym; \ + convolve_x_table[R][1] = &detail::convolve_x_radius_anti; \ + convolve_strided_table[R][0] = &detail::convolve_strided_radius_sym; \ + convolve_strided_table[R][1] = &detail::convolve_strided_radius_anti; \ + } + BIO_FILL(1) BIO_FILL(2) ... BIO_FILL(12) + #undef BIO_FILL +} +``` + +(Internally split each existing `template ` into +two non-templated-on-`Symmetric` aliases — `_sym` and `_anti` — so the +function-pointer types are concrete and the table is plain data.) + +Call `dispatch::init()` from `bind_filters()` in `src/bindings/filters.cxx` +(once, before any kernel binding can be invoked). Use a +`static std::once_flag` guard so re-imports don't double-initialise. + +### AVX2 kernel skeleton + +The X-pass kernel in `convolve_avx2.cxx` is essentially the scalar main +loop with explicit `__m256` registers: + +```cpp +namespace bioimage_cpp::filters::avx2 { + +template +void convolve_x_radius_sym( + const float* __restrict in, + float* __restrict out, + std::ptrdiff_t n_rows, + std::ptrdiff_t n_cols, + const float* __restrict h +) { + const std::ptrdiff_t prologue_end = std::min(R, n_cols); + const std::ptrdiff_t epilogue_start = std::max(prologue_end, n_cols - R); + + for (std::ptrdiff_t row = 0; row < n_rows; ++row) { + const float* __restrict in_row = in + row * n_cols; + float* __restrict out_row = out + row * n_cols; + + // --- border prologue: reuse scalar mirror code unchanged --- + scalar_border_sym(in_row, out_row, 0, prologue_end, n_cols, h); + + // --- main AVX2 loop --- + std::ptrdiff_t x = prologue_end; + const __m256 h0 = _mm256_set1_ps(h[0]); + for (; x + 8 <= epilogue_start; x += 8) { + __m256 acc = _mm256_mul_ps(_mm256_loadu_ps(in_row + x), h0); + for (int k = 1; k <= R; ++k) { + const __m256 hk = _mm256_set1_ps(h[k]); + const __m256 sum = _mm256_add_ps( + _mm256_loadu_ps(in_row + x + k), + _mm256_loadu_ps(in_row + x - k) + ); + acc = _mm256_fmadd_ps(hk, sum, acc); + } + _mm256_storeu_ps(out_row + x, acc); + } + // --- scalar tail (0..7 floats) --- + scalar_main_sym(in_row, out_row, x, epilogue_start, h); + + // --- border epilogue --- + scalar_border_sym(in_row, out_row, epilogue_start, n_cols, n_cols, h); + } +} + +template +void convolve_x_radius_anti(...) { /* same shape, _mm256_sub_ps instead of _add_ps */ } + +} +``` + +The strided kernel follows the same pattern but loops over `kStripBlock` +in steps of 8, using `__m256` for the accumulator strip. Crucially the +strip stays at 64 floats (`kStripBlock` is already a multiple of 8), so +no new tiling decision is needed. + +`scalar_border_sym` / `scalar_main_sym` are just the existing scalar +loop bodies hoisted into small inline helpers callable from both the +AVX2 and the scalar TU. **The mirror-boundary handling code must not be +duplicated between the two TUs** — that's where divergence bugs would +hide. Make the helpers `inline` in a shared header. + +### What stays byte-for-byte identical + +- Kernel-coefficient generation in `kernel.hxx`. +- Eigenvalue solvers in `eigenvalues.hxx`. +- Composite filters in `gaussian.hxx`. +- Binding layer in `src/bindings/filters.cxx`. +- Python wrapper in `src/bioimage_cpp/filters/_filters.py`. +- The public `convolve_axis_x` / `convolve_axis_strided` signatures in + `convolve.hxx`. +- The mirror-index function `detail::mirror_index` and the border + prologue/epilogue logic. + +If a Tier-2 change is touching any of these, stop and re-read the scope +section — it's almost certainly not what Tier 2 is for. + +### Expected speedup + +Based on the gap to `fastfilters` in this benchmark +(`bioimage_cpp / fastfilters` geomean = 2.00): + +- Simple filters (smoothing, derivative, gradient_magnitude, LoG): + realistic post-Tier-2 ratio **1.0 – 1.3×** fastfilters (essentially + tied to slightly behind). +- 3D Hessian / structure-tensor eigenvalues: realistic post-Tier-2 + ratio **1.4 – 1.7×** fastfilters. The remaining gap is in `acos`/`cos` + inside the 3×3 trig eigensolver, which intrinsics alone do not help + with. + +Do not expect to *match* fastfilters exactly without also vendoring a +vectorized math library and per-radius file-copy specialisation — that's +the next-tier-after-Tier-2 work, deliberately out of scope here. + +### Verification + +1. **Parity gate stays green**: + `python development/filters/check_parity.py` on a machine with AVX2 + support must still PASS at the same tolerances. If it doesn't, the + AVX2 kernel disagrees with the scalar kernel — that's the most + likely failure mode and is almost always an off-by-one in the + prologue / main / epilogue boundary handling. +2. **Scalar path stays green**: re-run with the AVX2 path forced off + (set the function-pointer table to the scalar entries unconditionally + in a debug build, or guard the dispatch decision with an environment + variable like `BIOIMAGE_FORCE_SCALAR=1`). The full pytest suite must + still pass; this catches scalar-only regressions introduced when + refactoring shared helpers. +3. **Benchmark**: + `python development/filters/benchmark.py` should show the + `bioimage_cpp / fastfilters` geomean drop from ~2.00 toward ~1.2. + Update this file with the new numbers. +4. **Pre-Haswell smoke**: the dispatch must take the scalar path on a + machine without AVX2. Easiest local check: temporarily make + `detect_avx2_fma()` return `false` and confirm correctness + + performance fall back to today's Tier-1 numbers. + +### Smallest first step + +Don't ship both kernels at once. The recommended sequence is: + +1. Land the dispatch scaffolding (`convolve_dispatch.hxx` / + `.cxx`, function-pointer tables, CMake wiring) **with both pointers + still pointing at the existing scalar kernels**. No behavior change. + Tests stay green. This isolates the build-system part of the work. +2. Add `convolve_x_radius_avx2` only. Re-run parity + benchmark. + Expect the simple filters to move; the Y/Z-bound filters + (gradient_magnitude, LoG, Hessian) move proportionally less. +3. Add `convolve_strided_radius_avx2`. Re-run parity + benchmark. + Expect the 3D filters to move significantly. + +If after step 2 the speedup is smaller than expected, stop and profile +before continuing — it usually means the autovectorizer was already +doing better than this section assumes, and the marginal value of +step 3 is lower than it appears here. + +### Watch-outs + +- **Boundary-mode divergence** between scalar and AVX2 paths is the + single most likely correctness bug. Share the prologue/epilogue + helpers via an `inline` header; don't copy-paste. +- **Unaligned loads only.** Use `_mm256_loadu_ps` / `_mm256_storeu_ps`, + not the aligned variants. The bench inputs are not guaranteed to be + 32-byte aligned, and on modern Intel/AMD the unaligned-load + performance penalty is essentially zero. Trying to force alignment in + the binding layer is more complexity than the win. +- **MSVC AVX2 detection.** `__builtin_cpu_supports` is GCC/Clang only. + Use raw `__cpuid` / `__cpuidex` + an `_xgetbv` check (the OS must + have saved the YMM state for AVX to be safe to use). There's example + code in `fastfilters/src/library/cpu_intel.c` if you need a + reference; do not vendor it, just write the small bit you need. +- **Don't introduce OpenMP, std::thread, or any threading primitive in + this work.** Threading is a separate follow-up that should layer on + top of the dispatch scheme via `parallel_for_chunks` (see + `include/bioimage_cpp/detail/threading.hxx`). Mixing the two changes + is asking for trouble. +- **Don't add AVX-512 "while we're here."** It is a separate trigger + decision with separate trade-offs (frequency throttling on older + Xeons, larger binary, marginal win on memory-bound separable FIR). + ## Known caveats reflected in the adapters - `fastfilters.gaussianDerivative` only accepts a uniform per-axis order; diff --git a/pyproject.toml b/pyproject.toml index ecbd508..b20b316 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -24,7 +24,7 @@ classifiers = [ ] [project.optional-dependencies] -test = ["pytest", "pooch"] +test = ["pytest", "pooch", "scipy"] data = ["pooch"] [tool.scikit-build]