From 92f5d352769463be09dbd457fb2b8548e898fb92 Mon Sep 17 00:00:00 2001 From: Constantin Pape Date: Sat, 15 Aug 2026 13:01:01 +0200 Subject: [PATCH] Expose structure tensor and implement eigenvector computation --- .github/workflows/wheels.yml | 2 +- MIGRATION_GUIDE.md | 15 + development/filters/PERFORMANCE_NOTES.md | 1 + .../filters/benchmark_structure_tensor.py | 337 ++++++++++++++++++ include/bioimage_cpp/filters/eigenvectors.hxx | 235 ++++++++++++ include/bioimage_cpp/filters/gaussian.hxx | 321 +++++++++++------ pyproject.toml | 2 +- src/bindings/filters.cxx | 229 ++++++++++++ src/bioimage_cpp/filters/__init__.py | 7 +- src/bioimage_cpp/filters/_filters.py | 100 +++++- tests/test_filters.py | 259 +++++++++++++- 11 files changed, 1386 insertions(+), 122 deletions(-) create mode 100644 development/filters/benchmark_structure_tensor.py create mode 100644 include/bioimage_cpp/filters/eigenvectors.hxx diff --git a/.github/workflows/wheels.yml b/.github/workflows/wheels.yml index 46cf80a..bab4a4c 100644 --- a/.github/workflows/wheels.yml +++ b/.github/workflows/wheels.yml @@ -27,7 +27,7 @@ jobs: CIBW_SKIP: "*-win32 *-manylinux_i686 *-musllinux_*" CIBW_ARCHS_MACOS: "x86_64 arm64" CIBW_ENVIRONMENT_MACOS: "MACOSX_DEPLOYMENT_TARGET=10.13" - CIBW_TEST_REQUIRES: "pytest pooch scipy" + CIBW_TEST_REQUIRES: "pytest pooch scipy scikit-image" CIBW_TEST_COMMAND: "python -c \"import bioimage_cpp; print(bioimage_cpp.__file__)\" && pytest -q {project}/tests" - name: Upload wheels diff --git a/MIGRATION_GUIDE.md b/MIGRATION_GUIDE.md index b3efa8b..1144eea 100644 --- a/MIGRATION_GUIDE.md +++ b/MIGRATION_GUIDE.md @@ -2259,6 +2259,10 @@ 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) + +tensor = bic.filters.structure_tensor(img, inner_sigma=1.0, outer_sigma=2.0) +# 3D component order: Jzz, Jzy, Jzx, Jyy, Jyx, Jxx +orientation = bic.filters.symmetric_eigenvector(tensor, index=2) ``` Name mapping: @@ -2270,6 +2274,7 @@ Name mapping: | `gaussianGradientMagnitude` | `gaussian_gradient_magnitude` | | `laplacianOfGaussian` | `laplacian_of_gaussian` | | `hessianOfGaussianEigenvalues` | `hessian_of_gaussian_eigenvalues` | +| `structureTensor` | `structure_tensor` | | `structureTensorEigenvalues` | `structure_tensor_eigenvalues` | Common parameters: @@ -2280,6 +2285,12 @@ Common parameters: per-axis sequence of ints in `{0, 1, 2}`. - `structure_tensor_eigenvalues` takes positional `inner_sigma` and `outer_sigma` (vigra calls them `innerScale` / `outerScale`). +- `structure_tensor` uses the same scales and returns a component-first array. + The 2D order is `(Jyy, Jyx, Jxx)`. The 3D order is + `(Jzz, Jzy, Jzx, Jyy, Jyx, Jxx)`. +- `symmetric_eigenvector(components, index, mask=None)` selects one + eigenvector. Indices use descending eigenvalue order. The optional boolean + mask sets unselected output vectors to zero. - `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 @@ -2300,6 +2311,10 @@ Important differences from vigra and fastfilters: - 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]`. +- `skimage.feature.structure_tensor` uses Sobel derivatives and one Gaussian + integration scale. `bioimage_cpp.filters.structure_tensor` instead uses + Gaussian derivatives at `inner_sigma` and smooths their products at + `outer_sigma`. - No IIR / recursive Gaussian, no `convolve` / `recursiveFilter2D`, no morphology, no nonlinear diffusion, and no non-local means in v1. Use `scipy.ndimage`, `skimage`, or the original diff --git a/development/filters/PERFORMANCE_NOTES.md b/development/filters/PERFORMANCE_NOTES.md index 42b469e..9eb2424 100644 --- a/development/filters/PERFORMANCE_NOTES.md +++ b/development/filters/PERFORMANCE_NOTES.md @@ -11,6 +11,7 @@ python development/filters/check_parity.py --force-scalar python development/filters/validate_eigenvalue_approximation.py python development/filters/benchmark_eigenvalues.py python development/filters/benchmark.py --repeats 5 +python development/filters/benchmark_structure_tensor.py --repeats 5 ``` The benchmark reports the median wall time across five interleaved calls after diff --git a/development/filters/benchmark_structure_tensor.py b/development/filters/benchmark_structure_tensor.py new file mode 100644 index 0000000..2ca6419 --- /dev/null +++ b/development/filters/benchmark_structure_tensor.py @@ -0,0 +1,337 @@ +"""Benchmark structure tensor components and selected eigenvectors. + +The SciPy baseline implements the same two-scale Gaussian-derivative tensor as +bioimage-cpp. The native skimage baseline pre-smooths the image, then uses +Sobel derivatives through ``skimage.feature.structure_tensor``. These two +baselines do not produce identical tensors. + +Run:: + + python development/filters/benchmark_structure_tensor.py --small + python development/filters/benchmark_structure_tensor.py --repeats 5 +""" + +from __future__ import annotations + +import argparse +import csv +import sys +from itertools import combinations_with_replacement + +import numpy as np +from scipy import ndimage +from skimage.feature import structure_tensor as skimage_structure_tensor + +from _bench_utils import ( + BenchConfig, + format_results_table, + load_2d, + load_3d, + time_interleaved, +) + + +LIBRARIES = ("bioimage_cpp", "scipy_numpy", "skimage_numpy") + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser( + description="Benchmark structure tensor components and selected eigenvectors." + ) + 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("--mask-threshold", type=float, default=0.5) + parser.add_argument("--repeats", type=int, default=5) + parser.add_argument("--small", action="store_true") + parser.add_argument("--no-2d", action="store_true") + parser.add_argument("--no-3d", action="store_true") + parser.add_argument("--csv") + return parser.parse_args() + + +def scipy_components(image: np.ndarray, cfg: BenchConfig) -> np.ndarray: + unit_orders = np.eye(image.ndim, dtype=int) + gradients = [ + ndimage.gaussian_filter( + image, + cfg.inner_sigma, + order=unit_orders[axis], + mode="mirror", + truncate=cfg.truncate, + ) + for axis in range(image.ndim) + ] + return np.stack( + [ + ndimage.gaussian_filter( + gradients[row] * gradients[column], + cfg.outer_sigma, + mode="mirror", + truncate=cfg.truncate, + ) + for row, column in combinations_with_replacement( + range(image.ndim), 2 + ) + ] + ) + + +def native_skimage_components(image: np.ndarray, cfg: BenchConfig) -> np.ndarray: + smoothed = ndimage.gaussian_filter( + image, + cfg.inner_sigma, + mode="mirror", + truncate=cfg.truncate, + ) + return np.ascontiguousarray( + skimage_structure_tensor( + smoothed, + sigma=cfg.outer_sigma, + mode="mirror", + order="rc", + ) + ) + + +def matrices_at_mask(components: np.ndarray, mask: np.ndarray) -> np.ndarray: + matrix_dimension = 2 if components.shape[0] == 3 else 3 + coordinates = np.nonzero(mask) + matrices = np.empty( + (coordinates[0].size, matrix_dimension, matrix_dimension), + dtype=components.dtype, + ) + for component, (row, column) in zip( + components, + combinations_with_replacement(range(matrix_dimension), 2), + strict=True, + ): + values = component[coordinates] + matrices[:, row, column] = values + matrices[:, column, row] = values + return matrices + + +def numpy_smallest_eigenvector( + components: np.ndarray, mask: np.ndarray +) -> np.ndarray: + matrix_dimension = 2 if components.shape[0] == 3 else 3 + matrices = matrices_at_mask(components, mask) + _, eigenvectors = np.linalg.eigh(matrices) + output = np.zeros(mask.shape + (matrix_dimension,), dtype=components.dtype) + output[mask] = eigenvectors[..., 0] + return output + + +def validate_exact_path( + image: np.ndarray, mask: np.ndarray, cfg: BenchConfig +) -> None: + from bioimage_cpp import filters as bf + + components = bf.structure_tensor( + image, + cfg.inner_sigma, + cfg.outer_sigma, + window_size=cfg.window_size, + ) + reference_components = scipy_components(image, cfg) + maximum_component_error = float( + np.max(np.abs(components - reference_components), initial=0.0) + ) + if maximum_component_error > 2e-3: + raise RuntimeError( + "structure tensor parity failed: maximum component error=" + f"{maximum_component_error:.3e}" + ) + + selected = bf.symmetric_eigenvector( + components, image.ndim - 1, mask=mask + ) + matrices = matrices_at_mask(components, mask) + eigenvalues, eigenvectors = np.linalg.eigh(matrices) + got = selected[mask] + expected = eigenvectors[..., 0] + dots = np.abs(np.sum(got * expected, axis=-1)) + scale = np.maximum( + np.linalg.norm(matrices, axis=(-2, -1)), + np.finfo(components.dtype).tiny, + ) + residual = np.linalg.norm( + np.einsum("...ij,...j->...i", matrices, got) + - eigenvalues[..., 0, None] * got, + axis=-1, + ) / scale + if np.max(residual, initial=0.0) > 2e-5: + raise RuntimeError( + "selected eigenvector residual failed: maximum residual=" + f"{np.max(residual):.3e}" + ) + + if matrices.shape[0] != 0: + gap = eigenvalues[..., 1] - eigenvalues[..., 0] + unique = gap > 1e-4 * scale + if np.any(unique) and np.min(dots[unique]) < 1.0 - 2e-4: + raise RuntimeError( + "selected eigenvector direction failed: minimum absolute dot=" + f"{np.min(dots[unique]):.6f}" + ) + print( + f" parity: max component error={maximum_component_error:.3e}, " + f"max eigenvector residual={np.max(residual, initial=0.0):.3e}" + ) + + +def summarize(results: dict[str, dict]) -> dict[str, dict]: + return { + library: {"median": result["median"], "min": result["min"]} + for library, result in results.items() + } + + +def benchmark_target( + label: str, + image: np.ndarray, + cfg: BenchConfig, + threshold: float, + repeats: int, +) -> list[dict]: + from bioimage_cpp import filters as bf + + mask = np.ascontiguousarray(image > threshold) + selected_index = image.ndim - 1 + print( + f"\n== {label}: shape={image.shape}, foreground={mask.mean():.1%} ==" + ) + validate_exact_path(image, mask, cfg) + + exact_components = scipy_components(image, cfg) + stages = { + "structure_tensor": { + "bioimage_cpp": lambda value: bf.structure_tensor( + value, + cfg.inner_sigma, + cfg.outer_sigma, + window_size=cfg.window_size, + ), + "scipy_numpy": lambda value: scipy_components(value, cfg), + "skimage_numpy": lambda value: native_skimage_components(value, cfg), + }, + "smallest_eigenvector": { + "bioimage_cpp": lambda _: bf.symmetric_eigenvector( + exact_components, selected_index, mask=mask + ), + "scipy_numpy": lambda _: numpy_smallest_eigenvector( + exact_components, mask + ), + }, + "end_to_end": { + "bioimage_cpp": lambda value: bf.symmetric_eigenvector( + bf.structure_tensor( + value, + cfg.inner_sigma, + cfg.outer_sigma, + window_size=cfg.window_size, + ), + selected_index, + mask=mask, + ), + "scipy_numpy": lambda value: numpy_smallest_eigenvector( + scipy_components(value, cfg), mask + ), + "skimage_numpy": lambda value: numpy_smallest_eigenvector( + native_skimage_components(value, cfg), mask + ), + }, + } + + rows = [] + for stage, callables in stages.items(): + results = time_interleaved(callables, image, repeats) + rows.append( + { + "filter": stage, + "dim": label, + "shape": str(tuple(image.shape)), + "results": summarize(results), + } + ) + return rows + + +def write_csv(path: str, rows: list[dict], repeats: int) -> None: + with open(path, "w", newline="") as file: + writer = csv.DictWriter( + file, + fieldnames=[ + "stage", "dim", "shape", "library", "median_s", "min_s", "repeats", + ], + ) + writer.writeheader() + for row in rows: + for library, result in row["results"].items(): + writer.writerow( + { + "stage": row["filter"], + "dim": row["dim"], + "shape": row["shape"], + "library": library, + "median_s": result["median"], + "min_s": result["min"], + "repeats": repeats, + } + ) + + +def main() -> int: + args = parse_args() + if args.repeats < 1: + raise ValueError("--repeats must be positive") + cfg = BenchConfig( + inner_sigma=args.inner_sigma, + outer_sigma=args.outer_sigma, + window_size=args.window_size, + ) + + targets = [] + if not args.no_2d: + targets.append( + ("2D", load_2d(crop=(128, 128) if args.small else None)) + ) + if not args.no_3d: + targets.append( + ("3D", load_3d(crop=(16, 64, 64) if args.small else None)) + ) + if not targets: + raise ValueError("enable at least one of the 2D and 3D benchmarks") + + print( + f"inner_sigma={cfg.inner_sigma}, outer_sigma={cfg.outer_sigma}, " + f"window_size={cfg.window_size}, threshold={args.mask_threshold}, " + f"repeats={args.repeats}" + ) + print( + "skimage_numpy uses Sobel derivatives and native skimage Gaussian " + "support; it is not a numerical-equivalence baseline." + ) + + rows = [] + for label, image in targets: + rows.extend( + benchmark_target( + label, + image, + cfg, + args.mask_threshold, + args.repeats, + ) + ) + print() + print(format_results_table(rows, libraries=LIBRARIES)) + if args.csv: + write_csv(args.csv, rows, args.repeats) + print(f"wrote {args.csv}") + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/include/bioimage_cpp/filters/eigenvectors.hxx b/include/bioimage_cpp/filters/eigenvectors.hxx new file mode 100644 index 0000000..ecd4630 --- /dev/null +++ b/include/bioimage_cpp/filters/eigenvectors.hxx @@ -0,0 +1,235 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +namespace bioimage_cpp::filters { + +namespace detail { + +template +inline void canonicalize_eigenvector(std::array &vector) { + std::size_t pivot = 0; + T largest = std::abs(vector[0]); + for (std::size_t axis = 1; axis < D; ++axis) { + const T magnitude = std::abs(vector[axis]); + if (magnitude > largest) { + largest = magnitude; + pivot = axis; + } + } + if (vector[pivot] < T{0}) { + for (auto &value : vector) { + value = -value; + } + } +} + +template +inline void jacobi_rotation( + std::array, D> &matrix, + std::array, D> &vectors, + const std::size_t p, + const std::size_t q +) { + const T off_diagonal = matrix[p][q]; + if (off_diagonal == T{0}) { + return; + } + + const T tau = (matrix[q][q] - matrix[p][p]) / (T{2} * off_diagonal); + const T sign = tau < T{0} ? T{-1} : T{1}; + const T tangent = sign / (std::abs(tau) + std::hypot(T{1}, tau)); + const T cosine = T{1} / std::sqrt(T{1} + tangent * tangent); + const T sine = tangent * cosine; + + const T diagonal_p = matrix[p][p]; + const T diagonal_q = matrix[q][q]; + matrix[p][p] = diagonal_p - tangent * off_diagonal; + matrix[q][q] = diagonal_q + tangent * off_diagonal; + matrix[p][q] = T{0}; + matrix[q][p] = T{0}; + + for (std::size_t axis = 0; axis < D; ++axis) { + if (axis == p || axis == q) { + continue; + } + const T value_p = matrix[axis][p]; + const T value_q = matrix[axis][q]; + const T rotated_p = cosine * value_p - sine * value_q; + const T rotated_q = sine * value_p + cosine * value_q; + matrix[axis][p] = rotated_p; + matrix[p][axis] = rotated_p; + matrix[axis][q] = rotated_q; + matrix[q][axis] = rotated_q; + } + + for (std::size_t axis = 0; axis < D; ++axis) { + const T value_p = vectors[axis][p]; + const T value_q = vectors[axis][q]; + vectors[axis][p] = cosine * value_p - sine * value_q; + vectors[axis][q] = sine * value_p + cosine * value_q; + } +} + +template +inline std::array symmetric_selected_eigenvector_descending( + std::array, D> matrix, + const std::size_t index +) { + T scale = T{0}; + for (const auto &row : matrix) { + for (const T value : row) { + scale = std::max(scale, std::abs(value)); + } + } + + if (scale == T{0}) { + std::array result{}; + result[index] = T{1}; + return result; + } + + if (scale < std::numeric_limits::min()) { + for (auto &row : matrix) { + for (auto &value : row) { + value /= scale; + } + } + } else { + const T inverse_scale = T{1} / scale; + for (auto &row : matrix) { + for (auto &value : row) { + value *= inverse_scale; + } + } + } + + std::array, D> vectors{}; + for (std::size_t axis = 0; axis < D; ++axis) { + vectors[axis][axis] = T{1}; + } + + constexpr std::size_t kMaxSweeps = 12; + const T tolerance = T{8} * std::numeric_limits::epsilon(); + for (std::size_t sweep = 0; sweep < kMaxSweeps; ++sweep) { + T off_diagonal_squared = T{0}; + for (std::size_t row = 0; row < D; ++row) { + for (std::size_t column = row + 1; column < D; ++column) { + off_diagonal_squared += + T{2} * matrix[row][column] * matrix[row][column]; + } + } + T diagonal_squared = T{0}; + for (std::size_t axis = 0; axis < D; ++axis) { + diagonal_squared += matrix[axis][axis] * matrix[axis][axis]; + } + if (std::sqrt(off_diagonal_squared) <= + tolerance * std::max(T{1}, std::sqrt(diagonal_squared))) { + break; + } + + for (std::size_t p = 0; p < D; ++p) { + for (std::size_t q = p + 1; q < D; ++q) { + jacobi_rotation(matrix, vectors, p, q); + } + } + } + + std::array order{}; + for (std::size_t axis = 0; axis < D; ++axis) { + order[axis] = axis; + } + std::stable_sort(order.begin(), order.end(), [&](const auto first, const auto second) { + return matrix[first][first] > matrix[second][second]; + }); + + std::array result{}; + T squared_norm = T{0}; + const std::size_t column = order[index]; + for (std::size_t axis = 0; axis < D; ++axis) { + result[axis] = vectors[axis][column]; + squared_norm += result[axis] * result[axis]; + } + const T inverse_norm = T{1} / std::sqrt(squared_norm); + for (auto &value : result) { + value *= inverse_norm; + } + canonicalize_eigenvector(result); + return result; +} + +} // namespace detail + +// Select one eigenvector from 2x2 symmetric matrices in descending eigenvalue +// order. Inputs use upper-triangle structure-of-arrays layout. +template +inline void ev2_symmetric_eigenvector_descending( + const T *__restrict a00, + const T *__restrict a01, + const T *__restrict a11, + const std::size_t index, + const std::uint8_t *__restrict mask, + T *__restrict out, + const std::ptrdiff_t n +) { + if (index >= 2) { + throw std::invalid_argument("eigenvector index must be in [0, 2)"); + } + for (std::ptrdiff_t i = 0; i < n; ++i) { + if (mask != nullptr && mask[i] == 0) { + out[2 * i] = T{0}; + out[2 * i + 1] = T{0}; + continue; + } + const auto vector = detail::symmetric_selected_eigenvector_descending( + {{{a00[i], a01[i]}, {a01[i], a11[i]}}}, index + ); + out[2 * i] = vector[0]; + out[2 * i + 1] = vector[1]; + } +} + +// Select one eigenvector from 3x3 symmetric matrices in descending eigenvalue +// order. Inputs use upper-triangle structure-of-arrays layout. +template +inline void ev3_symmetric_eigenvector_descending( + const T *__restrict a00, + const T *__restrict a01, + const T *__restrict a02, + const T *__restrict a11, + const T *__restrict a12, + const T *__restrict a22, + const std::size_t index, + const std::uint8_t *__restrict mask, + T *__restrict out, + const std::ptrdiff_t n +) { + if (index >= 3) { + throw std::invalid_argument("eigenvector index must be in [0, 3)"); + } + for (std::ptrdiff_t i = 0; i < n; ++i) { + if (mask != nullptr && mask[i] == 0) { + out[3 * i] = T{0}; + out[3 * i + 1] = T{0}; + out[3 * i + 2] = T{0}; + continue; + } + const auto vector = detail::symmetric_selected_eigenvector_descending( + {{{a00[i], a01[i], a02[i]}, + {a01[i], a11[i], a12[i]}, + {a02[i], a12[i], a22[i]}}}, + index + ); + out[3 * i] = vector[0]; + out[3 * i + 1] = vector[1]; + out[3 * i + 2] = vector[2]; + } +} + +} // namespace bioimage_cpp::filters diff --git a/include/bioimage_cpp/filters/gaussian.hxx b/include/bioimage_cpp/filters/gaussian.hxx index a23a43d..9946602 100644 --- a/include/bioimage_cpp/filters/gaussian.hxx +++ b/include/bioimage_cpp/filters/gaussian.hxx @@ -590,26 +590,26 @@ inline void hessian_of_gaussian_eigenvalues_3d( } // --------------------------------------------------------------------------- -// 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). +// Structure tensor. The leading output axis stores the upper triangle in +// NumPy axis order: (yy, yx, xx) or (zz, zy, zx, yy, yx, xx). // --------------------------------------------------------------------------- -inline void structure_tensor_eigenvalues_2d( +namespace detail { + +template +inline void structure_tensor_components_2d_profiled( 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::array &components, + const std::array &workspace, + const std::ptrdiff_t ny, + const std::ptrdiff_t nx, + const double sigma_inner_y, + const double sigma_inner_x, + const double sigma_outer_y, + const double sigma_outer_x, + const double window_ratio, + Profiler &profile ) { - BIOIMAGE_PROFILE_INIT(profile); - 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); @@ -617,70 +617,45 @@ inline void structure_tensor_eigenvalues_2d( const auto koy0 = gaussian_kernel(sigma_outer_y, 0, window_ratio); const auto kox0 = gaussian_kernel(sigma_outer_x, 0, window_ratio); - std::unique_ptr scratch; - { - BIOIMAGE_PROFILE_SCOPE(profile, "scratch_alloc"); - scratch = detail::allocate_scratch(n, 6); - } - float *work = detail::scratch_slot(scratch, n, 0); - float *gy = detail::scratch_slot(scratch, n, 1); - float *gx = detail::scratch_slot(scratch, n, 2); - float *syy_y = detail::scratch_slot(scratch, n, 3); - float *syx_y = detail::scratch_slot(scratch, n, 4); - float *sxx_y = detail::scratch_slot(scratch, n, 5); - - detail::gaussian_separable_2d_profiled(in, gy, work, ny, nx, kiy1, kix0, profile); - detail::gaussian_separable_2d_profiled(in, gx, work, ny, nx, kiy0, kix1, profile); + float *gy = components[0]; + float *gx = components[1]; + float *work = components[2]; + gaussian_separable_2d_profiled(in, gy, work, ny, nx, kiy1, kix0, profile); + gaussian_separable_2d_profiled(in, gx, work, ny, nx, kiy0, kix1, profile); { BIOIMAGE_PROFILE_SCOPE(profile, "outer_products"); const std::array gradients{gy, gx}; - const std::array components{syy_y, syx_y, sxx_y}; convolve_axis_strided_outer_products<2>( - gradients, components, 1, ny, nx, koy0 + gradients, workspace, 1, ny, nx, koy0 ); } { BIOIMAGE_PROFILE_SCOPE(profile, "axis_x"); - convolve_axis_x(syy_y, gy, ny, nx, kox0); - convolve_axis_x(syx_y, gx, ny, nx, kox0); - convolve_axis_x(sxx_y, work, ny, nx, kox0); + convolve_axis_x(workspace[0], components[0], ny, nx, kox0); + convolve_axis_x(workspace[1], components[1], ny, nx, kox0); + convolve_axis_x(workspace[2], components[2], ny, nx, kox0); } - - { - BIOIMAGE_PROFILE_SCOPE(profile, "eigenvalues"); - for (std::ptrdiff_t i = 0; i < n; ++i) { - const float a = gy[i]; - const float b = gx[i]; - const float c = work[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; - } - } - BIOIMAGE_PROFILE_REPORT(profile); } -inline void structure_tensor_eigenvalues_3d( +template +inline void structure_tensor_components_3d_profiled( 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::array &components, + const std::array &workspace, + const std::ptrdiff_t nz, + const std::ptrdiff_t ny, + const std::ptrdiff_t nx, + const double sigma_inner_z, + const double sigma_inner_y, + const double sigma_inner_x, + const double sigma_outer_z, + const double sigma_outer_y, + const double sigma_outer_x, + const double window_ratio, + Profiler &profile ) { - BIOIMAGE_PROFILE_INIT(profile); - 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); @@ -691,66 +666,208 @@ inline void structure_tensor_eigenvalues_3d( const auto koy0 = gaussian_kernel(sigma_outer_y, 0, window_ratio); const auto kox0 = gaussian_kernel(sigma_outer_x, 0, window_ratio); - std::unique_ptr scratch; - { - BIOIMAGE_PROFILE_SCOPE(profile, "scratch_alloc"); - scratch = detail::allocate_scratch(n, 10); - } - float *work = detail::scratch_slot(scratch, n, 0); - float *gz = detail::scratch_slot(scratch, n, 1); - float *gy = detail::scratch_slot(scratch, n, 2); - float *gx = detail::scratch_slot(scratch, n, 3); - float *szz = detail::scratch_slot(scratch, n, 4); - float *szy = detail::scratch_slot(scratch, n, 5); - float *szx = detail::scratch_slot(scratch, n, 6); - float *syy = detail::scratch_slot(scratch, n, 7); - float *syx = detail::scratch_slot(scratch, n, 8); - float *sxx = detail::scratch_slot(scratch, n, 9); - - detail::gaussian_first_axis_3d_profiled(in, gx, nz, ny, nx, kiz0, profile); - detail::gaussian_remaining_axes_3d_profiled( + float *work = workspace[0]; + float *gz = workspace[1]; + float *gy = workspace[2]; + float *gx = workspace[3]; + + gaussian_first_axis_3d_profiled(in, gx, nz, ny, nx, kiz0, profile); + gaussian_remaining_axes_3d_profiled( gx, gy, work, nz, ny, nx, kiy1, kix0, profile ); - detail::gaussian_remaining_axes_3d_profiled( + gaussian_remaining_axes_3d_profiled( gx, gx, work, nz, ny, nx, kiy0, kix1, profile ); - detail::gaussian_first_axis_3d_profiled(in, gz, nz, ny, nx, kiz1, profile); - detail::gaussian_remaining_axes_3d_profiled( + gaussian_first_axis_3d_profiled(in, gz, nz, ny, nx, kiz1, profile); + gaussian_remaining_axes_3d_profiled( gz, gz, work, nz, ny, nx, kiy0, kix0, profile ); { BIOIMAGE_PROFILE_SCOPE(profile, "outer_products"); const std::array gradients{gz, gy, gx}; - const std::array components{szz, szy, szx, syy, syx, sxx}; convolve_axis_strided_outer_products<3>( gradients, components, 1, nz, ny * nx, koz0 ); } - detail::gaussian_remaining_axes_3d_profiled( - szz, szz, work, nz, ny, nx, koy0, kox0, profile - ); - detail::gaussian_remaining_axes_3d_profiled( - szy, szy, work, nz, ny, nx, koy0, kox0, profile - ); - detail::gaussian_remaining_axes_3d_profiled( - szx, szx, work, nz, ny, nx, koy0, kox0, profile + for (float *component : components) { + gaussian_remaining_axes_3d_profiled( + component, component, work, nz, ny, nx, koy0, kox0, profile + ); + } +} + +} // namespace detail + +inline void structure_tensor_2d( + const float *in, + float *out, + const std::ptrdiff_t ny, + const std::ptrdiff_t nx, + const double sigma_inner_y, + const double sigma_inner_x, + const double sigma_outer_y, + const double sigma_outer_x, + const double window_ratio +) { + BIOIMAGE_PROFILE_INIT(profile); + const std::ptrdiff_t n = ny * nx; + const std::array components{out, out + n, out + 2 * n}; + std::unique_ptr scratch; + { + BIOIMAGE_PROFILE_SCOPE(profile, "scratch_alloc"); + scratch = detail::allocate_scratch(n, 3); + } + const std::array workspace{ + detail::scratch_slot(scratch, n, 0), + detail::scratch_slot(scratch, n, 1), + detail::scratch_slot(scratch, n, 2), + }; + detail::structure_tensor_components_2d_profiled( + in, components, workspace, ny, nx, + sigma_inner_y, sigma_inner_x, sigma_outer_y, sigma_outer_x, + window_ratio, profile ); - detail::gaussian_remaining_axes_3d_profiled( - syy, syy, work, nz, ny, nx, koy0, kox0, profile + BIOIMAGE_PROFILE_REPORT(profile); +} + +inline void structure_tensor_3d( + const float *in, + float *out, + const std::ptrdiff_t nz, + const std::ptrdiff_t ny, + const std::ptrdiff_t nx, + const double sigma_inner_z, + const double sigma_inner_y, + const double sigma_inner_x, + const double sigma_outer_z, + const double sigma_outer_y, + const double sigma_outer_x, + const double window_ratio +) { + BIOIMAGE_PROFILE_INIT(profile); + const std::ptrdiff_t n = nz * ny * nx; + const std::array components{ + out, out + n, out + 2 * n, out + 3 * n, out + 4 * n, out + 5 * n, + }; + std::unique_ptr scratch; + { + BIOIMAGE_PROFILE_SCOPE(profile, "scratch_alloc"); + scratch = detail::allocate_scratch(n, 4); + } + const std::array workspace{ + detail::scratch_slot(scratch, n, 0), + detail::scratch_slot(scratch, n, 1), + detail::scratch_slot(scratch, n, 2), + detail::scratch_slot(scratch, n, 3), + }; + detail::structure_tensor_components_3d_profiled( + in, components, workspace, nz, ny, nx, + sigma_inner_z, sigma_inner_y, sigma_inner_x, + sigma_outer_z, sigma_outer_y, sigma_outer_x, + window_ratio, profile ); - detail::gaussian_remaining_axes_3d_profiled( - syx, syx, work, nz, ny, nx, koy0, kox0, profile + BIOIMAGE_PROFILE_REPORT(profile); +} + +inline void structure_tensor_eigenvalues_2d( + const float *in, + float *out, + const std::ptrdiff_t ny, + const std::ptrdiff_t nx, + const double sigma_inner_y, + const double sigma_inner_x, + const double sigma_outer_y, + const double sigma_outer_x, + const double window_ratio +) { + BIOIMAGE_PROFILE_INIT(profile); + const std::ptrdiff_t n = ny * nx; + std::unique_ptr scratch; + { + BIOIMAGE_PROFILE_SCOPE(profile, "scratch_alloc"); + scratch = detail::allocate_scratch(n, 6); + } + const std::array components{ + detail::scratch_slot(scratch, n, 0), + detail::scratch_slot(scratch, n, 1), + detail::scratch_slot(scratch, n, 2), + }; + const std::array workspace{ + detail::scratch_slot(scratch, n, 3), + detail::scratch_slot(scratch, n, 4), + detail::scratch_slot(scratch, n, 5), + }; + detail::structure_tensor_components_2d_profiled( + in, components, workspace, ny, nx, + sigma_inner_y, sigma_inner_x, sigma_outer_y, sigma_outer_x, + window_ratio, profile ); - detail::gaussian_remaining_axes_3d_profiled( - sxx, sxx, work, nz, ny, nx, koy0, kox0, profile + + { + BIOIMAGE_PROFILE_SCOPE(profile, "eigenvalues"); + for (std::ptrdiff_t i = 0; i < n; ++i) { + const float half_trace = 0.5f * (components[0][i] + components[2][i]); + const float half_difference = + 0.5f * (components[0][i] - components[2][i]); + const float discriminant = std::sqrt( + half_difference * half_difference + components[1][i] * components[1][i] + ); + out[2 * i] = half_trace + discriminant; + out[2 * i + 1] = half_trace - discriminant; + } + } + BIOIMAGE_PROFILE_REPORT(profile); +} + +inline void structure_tensor_eigenvalues_3d( + const float *in, + float *out, + const std::ptrdiff_t nz, + const std::ptrdiff_t ny, + const std::ptrdiff_t nx, + const double sigma_inner_z, + const double sigma_inner_y, + const double sigma_inner_x, + const double sigma_outer_z, + const double sigma_outer_y, + const double sigma_outer_x, + const double window_ratio +) { + BIOIMAGE_PROFILE_INIT(profile); + const std::ptrdiff_t n = nz * ny * nx; + std::unique_ptr scratch; + { + BIOIMAGE_PROFILE_SCOPE(profile, "scratch_alloc"); + scratch = detail::allocate_scratch(n, 10); + } + const std::array workspace{ + detail::scratch_slot(scratch, n, 0), + detail::scratch_slot(scratch, n, 1), + detail::scratch_slot(scratch, n, 2), + detail::scratch_slot(scratch, n, 3), + }; + const std::array components{ + detail::scratch_slot(scratch, n, 4), + detail::scratch_slot(scratch, n, 5), + detail::scratch_slot(scratch, n, 6), + detail::scratch_slot(scratch, n, 7), + detail::scratch_slot(scratch, n, 8), + detail::scratch_slot(scratch, n, 9), + }; + detail::structure_tensor_components_3d_profiled( + in, components, workspace, nz, ny, nx, + sigma_inner_z, sigma_inner_y, sigma_inner_x, + sigma_outer_z, sigma_outer_y, sigma_outer_x, + window_ratio, profile ); { BIOIMAGE_PROFILE_SCOPE(profile, "eigenvalues"); ev3_symmetric_descending_interleaved( - szz, szy, szx, syy, syx, sxx, out, n + components[0], components[1], components[2], components[3], + components[4], components[5], out, n ); } BIOIMAGE_PROFILE_REPORT(profile); diff --git a/pyproject.toml b/pyproject.toml index 91dc4d4..07d8ec6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -34,7 +34,7 @@ Repository = "https://github.com/computational-cell-analytics/bioimage-cpp" Issues = "https://github.com/computational-cell-analytics/bioimage-cpp/issues" [project.optional-dependencies] -test = ["pytest", "pooch", "scipy"] +test = ["pytest", "pooch", "scipy", "scikit-image"] data = ["pooch"] [tool.scikit-build] diff --git a/src/bindings/filters.cxx b/src/bindings/filters.cxx index 0eefbf5..39a008d 100644 --- a/src/bindings/filters.cxx +++ b/src/bindings/filters.cxx @@ -1,15 +1,20 @@ #include "filters.hxx" #include "ndarray.hxx" +#include "bioimage_cpp/filters/eigenvectors.hxx" #include "bioimage_cpp/filters/gaussian.hxx" #include +#include #include #include +#include #include +#include #include #include +#include namespace nb = nanobind; @@ -18,6 +23,13 @@ namespace { using ConstImage = nb::ndarray; using Image = nb::ndarray; +using ConstMask = nb::ndarray; + +template +using ConstTypedArray = nb::ndarray; + +template +using TypedArray = nb::ndarray; void require_ndim(const ConstImage &image, int expected, const char *function) { if (static_cast(image.ndim()) != expected) { @@ -59,6 +71,12 @@ Image allocate_image(const std::size_t *shape, std::size_t ndim) { return detail::make_array(std::span(shape, ndim)); } +Image allocate_image_for_overwrite(const std::size_t *shape, std::size_t ndim) { + return detail::make_array_for_overwrite( + std::span(shape, ndim) + ); +} + // --------------------------------------------------------------------------- // gaussian_smoothing // --------------------------------------------------------------------------- @@ -391,6 +409,84 @@ Image hessian_of_gaussian_eigenvalues_3d( return out; } +// --------------------------------------------------------------------------- +// Structure tensor. Output has a leading upper-triangle component axis. +// --------------------------------------------------------------------------- + +Image structure_tensor_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_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] = {3, ny, nx}; + Image out = allocate_image_for_overwrite(shape, 3); + + const float *in_ptr = image.data(); + float *out_ptr = out.data(); + { + nb::gil_scoped_release release; + filters::structure_tensor_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_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_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] = {6, nz, ny, nx}; + Image out = allocate_image_for_overwrite(shape, 4); + + const float *in_ptr = image.data(); + float *out_ptr = out.data(); + { + nb::gil_scoped_release release; + filters::structure_tensor_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; +} + // --------------------------------------------------------------------------- // Structure-tensor eigenvalues. Output shape: input shape + (N,) trailing. // --------------------------------------------------------------------------- @@ -522,6 +618,110 @@ void ev3_symmetric_float32(ConstImage components, Image out) { } } +template +TypedArray symmetric_eigenvector( + ConstTypedArray components, + const std::size_t index, + std::optional mask +) { + const char *fn = "symmetric_eigenvector"; + if (components.ndim() < 1) { + throw std::invalid_argument( + std::string(fn) + ": components must have ndim >= 1" + ); + } + + const std::size_t component_count = components.shape(0); + std::size_t matrix_dimension = 0; + if (component_count == 3) { + matrix_dimension = 2; + } else if (component_count == 6) { + matrix_dimension = 3; + } else { + throw std::invalid_argument( + std::string(fn) + ": components.shape[0] must be 3 or 6, got " + + std::to_string(component_count) + ); + } + if (index >= matrix_dimension) { + throw std::invalid_argument( + std::string(fn) + ": index must be in [0, " + + std::to_string(matrix_dimension) + "), got " + std::to_string(index) + ); + } + + std::vector batch_shape(components.ndim() - 1); + for (std::size_t axis = 1; axis < components.ndim(); ++axis) { + batch_shape[axis - 1] = components.shape(axis); + } + const std::size_t n = detail::checked_array_size(batch_shape); + if (n > static_cast(std::numeric_limits::max())) { + throw std::invalid_argument(std::string(fn) + ": batch is too large"); + } + if (n != 0 && component_count > std::numeric_limits::max() / n) { + throw std::invalid_argument(std::string(fn) + ": component size overflows size_t"); + } + + const std::uint8_t *mask_ptr = nullptr; + if (mask.has_value()) { + if (mask->ndim() != batch_shape.size()) { + throw std::invalid_argument( + std::string(fn) + ": mask shape must match components batch shape" + ); + } + for (std::size_t axis = 0; axis < batch_shape.size(); ++axis) { + if (mask->shape(axis) != batch_shape[axis]) { + throw std::invalid_argument( + std::string(fn) + ": mask shape must match components batch shape" + ); + } + } + mask_ptr = mask->data(); + } + + const T *components_ptr = components.data(); + for (std::size_t i = 0; i < component_count * n; ++i) { + if (!std::isfinite(components_ptr[i])) { + throw std::invalid_argument( + std::string(fn) + ": components must contain only finite values" + ); + } + } + + std::vector output_shape = batch_shape; + output_shape.push_back(matrix_dimension); + auto output = detail::make_array_for_overwrite(output_shape); + T *output_ptr = output.data(); + { + nb::gil_scoped_release release; + if (matrix_dimension == 2) { + filters::ev2_symmetric_eigenvector_descending( + components_ptr, + components_ptr + n, + components_ptr + 2 * n, + index, + mask_ptr, + output_ptr, + static_cast(n) + ); + } else { + filters::ev3_symmetric_eigenvector_descending( + components_ptr, + components_ptr + n, + components_ptr + 2 * n, + components_ptr + 3 * n, + components_ptr + 4 * n, + components_ptr + 5 * n, + index, + mask_ptr, + output_ptr, + static_cast(n) + ); + } + } + return output; +} + } // namespace void bind_filters(nb::module_ &m) { @@ -538,6 +738,16 @@ void bind_filters(nb::module_ &m) { nb::arg("components"), nb::arg("out"), "Compute sorted eigenvalues for symmetric 3x3 float32 matrices." ); + m.def( + "_symmetric_eigenvector_float32", &symmetric_eigenvector, + nb::arg("components"), nb::arg("index"), nb::arg("mask") = nb::none(), + "Return one descending-order eigenvector from packed float32 matrices." + ); + m.def( + "_symmetric_eigenvector_float64", &symmetric_eigenvector, + nb::arg("components"), nb::arg("index"), nb::arg("mask") = nb::none(), + "Return one descending-order eigenvector from packed float64 matrices." + ); m.def( "_gaussian_smoothing_2d_float32", &gaussian_smoothing_2d, nb::arg("image"), nb::arg("sigma_y"), nb::arg("sigma_x"), @@ -607,6 +817,25 @@ void bind_filters(nb::module_ &m) { "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_2d_float32", &structure_tensor_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_size") = 0.0, + "Structure tensor of a 2D float32 image. " + "Output shape: (3, ny, nx), ordered as (Jyy, Jyx, Jxx)." + ); + m.def( + "_structure_tensor_3d_float32", &structure_tensor_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_size") = 0.0, + "Structure tensor of a 3D float32 image. " + "Output shape: (6, nz, ny, nx), ordered as " + "(Jzz, Jzy, Jzx, Jyy, Jyx, Jxx)." + ); m.def( "_structure_tensor_eigenvalues_2d_float32", &structure_tensor_eigenvalues_2d, nb::arg("image"), diff --git a/src/bioimage_cpp/filters/__init__.py b/src/bioimage_cpp/filters/__init__.py index 716d3fe..454e4c2 100644 --- a/src/bioimage_cpp/filters/__init__.py +++ b/src/bioimage_cpp/filters/__init__.py @@ -1,5 +1,4 @@ -"""Image filters: separable Gaussian-family derivatives, gradient magnitude, -Laplacian of Gaussian, Hessian and structure-tensor eigenvalues.""" +"""Gaussian-family image filters and symmetric tensor operations.""" from ._filters import ( gaussian_derivative, @@ -7,7 +6,9 @@ gaussian_smoothing, hessian_of_gaussian_eigenvalues, laplacian_of_gaussian, + structure_tensor, structure_tensor_eigenvalues, + symmetric_eigenvector, ) __all__ = [ @@ -16,5 +17,7 @@ "gaussian_gradient_magnitude", "laplacian_of_gaussian", "hessian_of_gaussian_eigenvalues", + "structure_tensor", "structure_tensor_eigenvalues", + "symmetric_eigenvector", ] diff --git a/src/bioimage_cpp/filters/_filters.py b/src/bioimage_cpp/filters/_filters.py index 813660b..564adc4 100644 --- a/src/bioimage_cpp/filters/_filters.py +++ b/src/bioimage_cpp/filters/_filters.py @@ -12,8 +12,9 @@ ``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. +* Eigenvalue functions return a trailing axis of size ``image.ndim`` in + descending order. +* ``structure_tensor`` stores upper-triangle components on a leading axis. """ from __future__ import annotations @@ -229,6 +230,41 @@ def hessian_of_gaussian_eigenvalues( return _finalise(result, out_dtype) +def structure_tensor( + image: np.ndarray, + inner_sigma: float | Sequence[float], + outer_sigma: float | Sequence[float], + *, + window_size: float = 0.0, +) -> np.ndarray: + """Compute the upper-triangle components of the structure tensor. + + The result has a leading component axis. The 2D order is + ``(Jyy, Jyx, Jxx)``. The 3D order is + ``(Jzz, Jzy, Jzx, Jyy, Jyx, Jxx)``. + """ + function = "structure_tensor" + 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) + window = _normalize_window(window_size, function) + if prepared.ndim == 2: + result = _core._structure_tensor_2d_float32( + prepared, + inner[0], inner[1], + outer[0], outer[1], + window, + ) + else: + result = _core._structure_tensor_3d_float32( + prepared, + inner[0], inner[1], inner[2], + outer[0], outer[1], outer[2], + window, + ) + return _finalise(result, out_dtype) + + def structure_tensor_eigenvalues( image: np.ndarray, inner_sigma: float | Sequence[float], @@ -261,3 +297,63 @@ def structure_tensor_eigenvalues( window, ) return _finalise(result, out_dtype) + + +def symmetric_eigenvector( + components: np.ndarray, + index: int, + *, + mask: np.ndarray | None = None, +) -> np.ndarray: + """Return one eigenvector from packed symmetric matrices. + + ``components`` must have shape ``(3, *batch_shape)`` for 2D matrices or + ``(6, *batch_shape)`` for 3D matrices. The component orders are + ``(A00, A01, A11)`` and ``(A00, A01, A02, A11, A12, A22)``. Eigenvector + indices use descending eigenvalue order. + + The result has shape ``batch_shape + (matrix_dimension,)``. Entries + outside ``mask`` are zero. A repeated eigenspace has no unique direction; + the function returns a deterministic unit vector in that eigenspace. + """ + function = "symmetric_eigenvector" + array = np.asarray(components) + if array.ndim < 1: + raise ValueError(f"{function}: components must have ndim >= 1") + if array.shape[0] == 3: + matrix_dimension = 2 + elif array.shape[0] == 6: + matrix_dimension = 3 + else: + raise ValueError( + f"{function}: components.shape[0] must be 3 or 6, " + f"got {array.shape[0]}" + ) + selected = strict_index( + index, "index", minimum=0, maximum=matrix_dimension - 1 + ) + if array.dtype not in (np.dtype(np.float32), np.dtype(np.float64)): + raise TypeError( + f"{function}: components dtype must be float32 or float64, " + f"got dtype={array.dtype}" + ) + prepared = np.ascontiguousarray(array) + + mask_arg = None + if mask is not None: + mask_array = np.asarray(mask) + if mask_array.shape != prepared.shape[1:]: + raise ValueError( + f"{function}: mask shape must match components batch shape, " + f"got mask shape={mask_array.shape}, " + f"batch shape={prepared.shape[1:]}" + ) + if mask_array.dtype != np.dtype(bool): + raise TypeError( + f"{function}: mask must have dtype bool, got dtype={mask_array.dtype}" + ) + mask_arg = np.ascontiguousarray(mask_array.view(np.uint8)) + + if prepared.dtype == np.dtype(np.float32): + return _core._symmetric_eigenvector_float32(prepared, selected, mask_arg) + return _core._symmetric_eigenvector_float64(prepared, selected, mask_arg) diff --git a/tests/test_filters.py b/tests/test_filters.py index fe08a56..0013862 100644 --- a/tests/test_filters.py +++ b/tests/test_filters.py @@ -7,10 +7,12 @@ """ from concurrent.futures import ThreadPoolExecutor +from itertools import combinations_with_replacement import numpy as np import pytest from scipy import ndimage +from skimage.feature import structure_tensor as skimage_structure_tensor from bioimage_cpp import _core import bioimage_cpp.filters as bf @@ -179,12 +181,20 @@ def test_hessian_eigenvalues_descending_order_3d(): # Structure tensor eigenvalues # --------------------------------------------------------------------------- +def _structure_tensor_reference_2d(img, inner, outer, *, truncate=None): + kwargs = {"mode": "mirror"} + if truncate is not None: + kwargs["truncate"] = truncate + gy = ndimage.gaussian_filter(img, inner, order=[1, 0], **kwargs) + gx = ndimage.gaussian_filter(img, inner, order=[0, 1], **kwargs) + syy = ndimage.gaussian_filter(gy * gy, outer, **kwargs) + syx = ndimage.gaussian_filter(gy * gx, outer, **kwargs) + sxx = ndimage.gaussian_filter(gx * gx, outer, **kwargs) + return np.stack([syy, syx, sxx]) + + 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") + syy, syx, sxx = _structure_tensor_reference_2d(img, inner, outer) mat = np.stack( [np.stack([syy, syx], axis=-1), np.stack([syx, sxx], axis=-1)], axis=-2, @@ -200,16 +210,26 @@ def test_structure_tensor_eigenvalues_2d_matches_reference(): np.testing.assert_allclose(got, ref, atol=2e-3) +def _structure_tensor_reference_3d(vol, inner, outer, *, truncate=None): + kwargs = {"mode": "mirror"} + if truncate is not None: + kwargs["truncate"] = truncate + gz = ndimage.gaussian_filter(vol, inner, order=[1, 0, 0], **kwargs) + gy = ndimage.gaussian_filter(vol, inner, order=[0, 1, 0], **kwargs) + gx = ndimage.gaussian_filter(vol, inner, order=[0, 0, 1], **kwargs) + szz = ndimage.gaussian_filter(gz * gz, outer, **kwargs) + szy = ndimage.gaussian_filter(gz * gy, outer, **kwargs) + szx = ndimage.gaussian_filter(gz * gx, outer, **kwargs) + syy = ndimage.gaussian_filter(gy * gy, outer, **kwargs) + syx = ndimage.gaussian_filter(gy * gx, outer, **kwargs) + sxx = ndimage.gaussian_filter(gx * gx, outer, **kwargs) + return np.stack([szz, szy, szx, syy, syx, sxx]) + + def _structure_tensor_eigenvalues_reference_3d(vol, inner, outer): - gz = ndimage.gaussian_filter(vol, inner, order=[1, 0, 0], mode="mirror") - gy = ndimage.gaussian_filter(vol, inner, order=[0, 1, 0], mode="mirror") - gx = ndimage.gaussian_filter(vol, inner, order=[0, 0, 1], mode="mirror") - szz = ndimage.gaussian_filter(gz * gz, outer, mode="mirror") - szy = ndimage.gaussian_filter(gz * gy, outer, mode="mirror") - szx = ndimage.gaussian_filter(gz * gx, outer, 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") + szz, szy, szx, syy, syx, sxx = _structure_tensor_reference_3d( + vol, inner, outer + ) mat = np.stack([ np.stack([szz, szy, szx], axis=-1), np.stack([szy, syy, syx], axis=-1), @@ -370,6 +390,191 @@ def test_structure_tensor_eigenvalues_3d_matches_reference(): assert np.all(got >= -1e-6) +@pytest.mark.parametrize( + ("shape", "inner", "outer", "reference"), + [ + ((23, 31), (0.8, 1.2), (1.3, 1.8), _structure_tensor_reference_2d), + ((7, 11, 13), (0.8, 1.0, 1.2), (1.3, 1.5, 1.7), _structure_tensor_reference_3d), + ], +) +def test_structure_tensor_components_match_scipy(shape, inner, outer, reference): + image = _random_image(shape) + got = bf.structure_tensor(image, inner, outer, window_size=3.0) + expected = reference(image, inner, outer, truncate=3.0) + component_count = image.ndim * (image.ndim + 1) // 2 + assert got.shape == (component_count,) + image.shape + assert got.dtype == np.float32 + np.testing.assert_allclose(got, expected, atol=2e-3) + + +def _matrices_from_components(components): + matrix_dimension = 2 if components.shape[0] == 3 else 3 + matrices = np.zeros( + components.shape[1:] + (matrix_dimension, matrix_dimension), + dtype=components.dtype, + ) + for component, (row, column) in zip( + components, + combinations_with_replacement(range(matrix_dimension), 2), + strict=True, + ): + matrices[..., row, column] = component + matrices[..., column, row] = component + return matrices + + +def _components_from_matrices(matrices): + matrix_dimension = matrices.shape[-1] + return np.ascontiguousarray( + np.stack( + [ + matrices[..., row, column] + for row, column in combinations_with_replacement( + range(matrix_dimension), 2 + ) + ] + ) + ) + + +def _assert_selected_eigenvector(matrices, got, index, *, rtol): + matrix_dimension = matrices.shape[-1] + eigenvalues, eigenvectors = np.linalg.eigh(matrices) + ascending_index = matrix_dimension - 1 - index + selected_value = eigenvalues[..., ascending_index] + expected = eigenvectors[..., :, ascending_index] + matrices64 = matrices.astype(np.float64) + got64 = got.astype(np.float64) + dots = np.abs(np.sum(got64 * expected.astype(np.float64), axis=-1)) + scale = np.maximum( + np.linalg.norm(matrices64, axis=(-2, -1)), + np.finfo(matrices.dtype).tiny, + ) + residual = np.linalg.norm( + np.einsum("...ij,...j->...i", matrices64, got64) + - selected_value.astype(np.float64)[..., None] * got64, + axis=-1, + ) / scale + np.testing.assert_allclose(np.linalg.norm(got64, axis=-1), 1.0, rtol=rtol) + assert np.max(residual, initial=0.0) < rtol + assert np.min(dots, initial=1.0) > 1.0 - 10.0 * rtol + + +@pytest.mark.parametrize("matrix_dimension", [2, 3]) +@pytest.mark.parametrize("dtype", [np.float32, np.float64]) +def test_symmetric_eigenvector_matches_numpy(matrix_dimension, dtype): + rng = np.random.default_rng(314 + matrix_dimension) + matrices = rng.normal(size=(5, 7, matrix_dimension, matrix_dimension)).astype(dtype) + matrices += matrices.swapaxes(-1, -2) + components = _components_from_matrices(matrices) + tolerance = 2e-5 if dtype == np.float32 else 1e-12 + for index in range(matrix_dimension): + got = bf.symmetric_eigenvector(components, index) + assert got.shape == matrices.shape[:-1] + assert got.dtype == dtype + _assert_selected_eigenvector(matrices, got, index, rtol=tolerance) + + +@pytest.mark.parametrize("shape", [(19, 23), (5, 7, 9)]) +def test_symmetric_eigenvector_accepts_skimage_components(shape): + image = _random_image(shape) + components = np.ascontiguousarray( + skimage_structure_tensor(image, sigma=1.2, mode="mirror", order="rc") + ) + matrices = _matrices_from_components(components) + for index in range(image.ndim): + got = bf.symmetric_eigenvector(components, index) + _assert_selected_eigenvector(matrices, got, index, rtol=2e-5) + + +def test_symmetric_eigenvector_mask_and_general_batch_shape(): + rng = np.random.default_rng(41) + matrices = rng.normal(size=(4, 5, 3, 3)).astype(np.float32) + matrices += matrices.swapaxes(-1, -2) + components = _components_from_matrices(matrices) + full = bf.symmetric_eigenvector(components, 2) + storage = rng.random((8, 10)) > 0.5 + mask = storage[::2, ::2] + assert not mask.flags.c_contiguous + got = bf.symmetric_eigenvector(components, 2, mask=mask) + np.testing.assert_array_equal(got[mask], full[mask]) + np.testing.assert_array_equal(got[~mask], 0.0) + + single = bf.symmetric_eigenvector(components[:, 0, 0], 1) + assert single.shape == (3,) + np.testing.assert_allclose(np.linalg.norm(single), 1.0, atol=2e-6) + + +@pytest.mark.parametrize("matrix_dimension", [2, 3]) +def test_symmetric_eigenvector_repeated_eigenspaces_are_valid(matrix_dimension): + rng = np.random.default_rng(73 + matrix_dimension) + matrices = [ + np.zeros((matrix_dimension, matrix_dimension)), + np.eye(matrix_dimension) * 4.0, + ] + q, _ = np.linalg.qr(rng.normal(size=(matrix_dimension, matrix_dimension))) + spectrum = np.ones(matrix_dimension) + spectrum[0] = 3.0 + matrices.append(q @ np.diag(spectrum) @ q.T) + matrices = np.asarray(matrices, dtype=np.float64) + components = _components_from_matrices(matrices) + + for index in range(matrix_dimension): + first = bf.symmetric_eigenvector(components, index) + second = bf.symmetric_eigenvector(components, index) + np.testing.assert_array_equal(first, second) + np.testing.assert_allclose(np.linalg.norm(first, axis=-1), 1.0, atol=1e-12) + eigenvalues = np.linalg.eigvalsh(matrices)[..., ::-1] + residual = np.linalg.norm( + np.einsum("...ij,...j->...i", matrices, first) + - eigenvalues[..., index, None] * first, + axis=-1, + ) + assert np.max(residual) < 1e-12 + pivot = np.argmax(np.abs(first), axis=-1) + assert np.all(np.take_along_axis(first, pivot[:, None], axis=-1) >= 0.0) + + +def test_symmetric_eigenvector_handles_float32_scales_and_near_repeated_roots(): + rng = np.random.default_rng(97) + matrices = [] + for scale in (np.finfo(np.float32).tiny, 1e-20, 1e-8, 1.0, 1e8, 1e20): + q, _ = np.linalg.qr(rng.normal(size=(3, 3))) + spectrum = np.asarray( + [3.0 * scale, scale * (1.0 + 2e-5), scale], dtype=np.float32 + ) + matrices.append((q @ np.diag(spectrum) @ q.T).astype(np.float32)) + matrices = np.stack(matrices) + components = _components_from_matrices(matrices) + for index in range(3): + got = bf.symmetric_eigenvector(components, index) + _assert_selected_eigenvector(matrices, got, index, rtol=3e-4) + + +def test_symmetric_eigenvector_validates_arguments(): + components = np.zeros((6, 4, 5), dtype=np.float32) + with pytest.raises(ValueError, match=r"shape\[0\]"): + bf.symmetric_eigenvector(np.zeros((5, 4), dtype=np.float32), 0) + with pytest.raises(ValueError, match="index"): + bf.symmetric_eigenvector(components, 3) + with pytest.raises(TypeError, match="integer"): + bf.symmetric_eigenvector(components, True) + with pytest.raises(TypeError, match="float32 or float64"): + bf.symmetric_eigenvector(components.astype(np.int32), 0) + with pytest.raises(ValueError, match="finite"): + invalid = components.copy() + invalid[0, 0, 0] = np.nan + bf.symmetric_eigenvector(invalid, 0) + with pytest.raises(ValueError, match="mask shape"): + bf.symmetric_eigenvector( + components, 0, mask=np.ones((4, 4), dtype=bool) + ) + with pytest.raises(TypeError, match="dtype bool"): + bf.symmetric_eigenvector( + components, 0, mask=np.ones((4, 5), dtype=np.uint8) + ) + + @pytest.mark.parametrize("shape", [(1, 5, 7), (2, 3, 4)]) def test_eigenvalue_filters_support_short_axes(shape): vol = _random_image(shape) @@ -402,6 +607,7 @@ def test_forced_scalar_matches_automatic_backend(monkeypatch): lambda: bf.gaussian_gradient_magnitude(vol, 1.5), lambda: bf.laplacian_of_gaussian(vol, 1.5), lambda: bf.hessian_of_gaussian_eigenvalues(vol, 1.5), + lambda: bf.structure_tensor(vol, 1.0, 2.0), lambda: bf.structure_tensor_eigenvalues(vol, 1.0, 2.0), ] automatic = [function() for function in functions] @@ -454,11 +660,29 @@ def test_concurrent_filter_calls_are_deterministic(): np.testing.assert_array_equal(result, expected) +def test_concurrent_symmetric_eigenvector_calls_are_deterministic(): + components = bf.structure_tensor(_random_image((7, 9, 11)), 1.0, 1.5) + expected = bf.symmetric_eigenvector(components, 2) + with ThreadPoolExecutor(max_workers=4) as executor: + results = list( + executor.map( + lambda _: bf.symmetric_eigenvector(components, 2), + range(8), + ) + ) + for result in results: + np.testing.assert_array_equal(result, expected) + + @pytest.mark.parametrize("shape", [(0, 4), (2, 0, 3)]) def test_filters_support_empty_inputs(shape): image = np.empty(shape, dtype=np.float32) assert bf.gaussian_smoothing(image, 1.0).shape == shape assert bf.hessian_of_gaussian_eigenvalues(image, 1.0).shape == shape + (len(shape),) + component_count = len(shape) * (len(shape) + 1) // 2 + components = bf.structure_tensor(image, 1.0, 2.0) + assert components.shape == (component_count,) + shape + assert bf.symmetric_eigenvector(components, 0).shape == shape + (len(shape),) assert ( bf.structure_tensor_eigenvalues(image, 1.0, 2.0).shape == shape + (len(shape),) @@ -476,6 +700,10 @@ def test_float64_input_returns_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) + tensor = bf.structure_tensor(img, 1.0, 2.0) + assert tensor.dtype == np.float64 + assert bf.symmetric_eigenvector(tensor, 1).dtype == np.float64 + def test_uint8_input_returns_float32(): rng = np.random.RandomState(0) @@ -502,6 +730,9 @@ def test_non_contiguous_input_is_handled(): 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) + tensor = bf.structure_tensor(sliced, 1.0, 2.0) + tensor_ref = _structure_tensor_reference_2d(sliced, 1.0, 2.0) + np.testing.assert_allclose(tensor, tensor_ref, atol=2e-3) # ---------------------------------------------------------------------------