From 2d3f14aa812f15b02342b0f4b416e44df4560cbd Mon Sep 17 00:00:00 2001 From: Constantin Pape Date: Tue, 4 Aug 2026 18:06:07 +0200 Subject: [PATCH 1/3] Implement ball invalidation and fix-branching off for TEASAR --- MIGRATION_GUIDE.md | 32 ++- development/skeleton/benchmark_teasar.py | 2 + .../skeleton/detail/compact_grid_dijkstra.hxx | 57 ++++ .../skeleton/detail/invalidation.hxx | 142 ++++++++++ include/bioimage_cpp/skeleton/teasar.hxx | 253 +++++++++++++----- src/bindings/skeleton.cxx | 29 +- src/bindings/skeleton_distributed.cxx | 20 +- src/bioimage_cpp/skeleton/__init__.py | 39 ++- src/bioimage_cpp/skeleton/distributed.py | 15 +- tests/skeleton/test_distributed.py | 46 +++- tests/skeleton/test_teasar.py | 100 ++++++- tests/skeleton/test_teasar_labels.py | 26 +- 12 files changed, 661 insertions(+), 100 deletions(-) create mode 100644 include/bioimage_cpp/skeleton/detail/invalidation.hxx diff --git a/MIGRATION_GUIDE.md b/MIGRATION_GUIDE.md index 5092fd1..b3efa8b 100644 --- a/MIGRATION_GUIDE.md +++ b/MIGRATION_GUIDE.md @@ -2635,6 +2635,8 @@ vertices, edges, radii = bic.skeleton.teasar( constant=0, pdrf_scale=100000, pdrf_exponent=4, + invalidation="ball", + fix_branching=True, number_of_threads=4, ) @@ -2646,6 +2648,8 @@ skeletons = bic.skeleton.teasar_labels( constant=0, pdrf_scale=100000, pdrf_exponent=4, + invalidation="ball", + fix_branching=True, number_of_threads=4, ) # {original_label: (vertices, edges, radii), ...} @@ -2689,14 +2693,20 @@ Important differences and current scope: distance-to-boundary field, a deterministic two-sweep root, a physical Dijkstra distance-from-root field, penalized repeated Dijkstra paths, and rolling invalidation with radius `scale * radius + constant`. -- This correctness-oriented implementation uses a physical axis-aligned - invalidation cube. The ordinary in-core entry points do not automatically - perform block stitching or accept manual targets; the dedicated block APIs - below provide mandatory interface targets and exact consolidation. Soma - handling, hole filling, component dust filtering, cross-section metadata, - and kimimaro's other postprocessing heuristics are not implemented. These - differences can change branch positions and vertex counts, so output is not - expected to be vertex-for-vertex identical to kimimaro. +- `invalidation="cube"` uses the existing physical axis-aligned invalidation + box and remains the default. `invalidation="ball"` uses an open physical + ball. It expands through active foreground with 26-connectivity, as in + kimimaro's component-aware ball invalidation. +- `fix_branching=True` remains the default. It routes each new path to the + existing zero-cost skeleton. Set `fix_branching=False` to build one parental + field from the root and extract all paths from that fixed field. +- The ordinary in-core entry points do not automatically perform block + stitching or accept manual targets. The dedicated block APIs below provide + mandatory interface targets and exact consolidation. Soma handling, hole + filling, component dust filtering, cross-section metadata, and kimimaro's + other postprocessing heuristics are not implemented. These differences can + change branch positions and vertex counts. Output is not expected to be + vertex-for-vertex identical to kimimaro. - The C++ core remains dependency-free. Component discovery uses x-runs and a union-find rather than dense component-label images. `number_of_threads=1` is the default and `0` uses hardware concurrency; one shared budget covers @@ -2746,6 +2756,8 @@ left_fragment = dist.block_teasar( origin=left_origin, required_targets=targets, spacing=spacing, + invalidation="cube", + fix_branching=True, ) # Repeat independently for every processing block. Neighboring calls select @@ -2779,6 +2791,10 @@ skeletons are connected through their anchors but are not expected to be vertex-for-vertex identical to whole-volume TEASAR because each block still has less path-selection context. +`block_teasar` and `block_teasar_labels` accept the same `invalidation` and +`fix_branching` options as the in-core entry points. Both defaults preserve the +previous block behavior. + Correctness tests are under `tests/skeleton/test_teasar.py`, `tests/skeleton/test_teasar_labels.py`, and `tests/skeleton/test_distributed.py`. The independent diff --git a/development/skeleton/benchmark_teasar.py b/development/skeleton/benchmark_teasar.py index 192282b..232f479 100644 --- a/development/skeleton/benchmark_teasar.py +++ b/development/skeleton/benchmark_teasar.py @@ -288,6 +288,8 @@ def bic_backend_call(mask, spacing, parameters, backend, number_of_threads=1): parameters["constant"], parameters["pdrf_scale"], parameters["pdrf_exponent"], + False, + True, backend, number_of_threads, ) diff --git a/include/bioimage_cpp/skeleton/detail/compact_grid_dijkstra.hxx b/include/bioimage_cpp/skeleton/detail/compact_grid_dijkstra.hxx index c156f10..7806b05 100644 --- a/include/bioimage_cpp/skeleton/detail/compact_grid_dijkstra.hxx +++ b/include/bioimage_cpp/skeleton/detail/compact_grid_dijkstra.hxx @@ -340,6 +340,63 @@ inline void compact_physical_distance_field( } } +template +inline void compact_node_cost_parental_field( + const CompactGridDomain &domain, + const std::uint32_t source, + const std::vector &costs, + CompactDijkstraWorkspace &workspace, + std::vector &predecessors, + CompactDijkstraStats *stats = nullptr +) { + const auto n = domain.size(); + if (source >= n || costs.size() != n) { + throw std::invalid_argument("invalid compact node-cost parental-field inputs"); + } + if constexpr (Adjacency == CompactAdjacency::Csr) { + if (!domain.has_csr()) { + throw std::invalid_argument("compact CSR adjacency is not available"); + } + } else if (!domain.has_full_lookup()) { + throw std::invalid_argument("compact full-index lookup is not available"); + } + + predecessors.resize(n); + workspace.state.assign(n, 0); + workspace.heap.clear(); + if (stats != nullptr) { + stats->reset(); + } + workspace.state[source] = kCompactDiscovered; + predecessors[source] = source; + compact_heap_push(workspace, {Distance{0}, source}, stats); + + while (!workspace.heap.empty()) { + const auto entry = compact_heap_pop(workspace, stats); + const auto node = entry.node; + if ((workspace.state[node] & kCompactSettled) != 0) { + continue; + } + workspace.state[node] |= kCompactSettled; + for_each_compact_neighbor( + domain, node, + [&](const std::uint32_t target, const double) { + if ((workspace.state[target] & + (kCompactDiscovered | kCompactSettled)) != 0) { + return; + } + workspace.state[target] |= kCompactDiscovered; + predecessors[target] = node; + compact_heap_push( + workspace, + {static_cast(entry.distance + costs[target]), target}, + stats + ); + } + ); + } +} + template inline void compact_node_cost_path( const CompactGridDomain &domain, diff --git a/include/bioimage_cpp/skeleton/detail/invalidation.hxx b/include/bioimage_cpp/skeleton/detail/invalidation.hxx new file mode 100644 index 0000000..14fe0b8 --- /dev/null +++ b/include/bioimage_cpp/skeleton/detail/invalidation.hxx @@ -0,0 +1,142 @@ +#pragma once + +#include "bioimage_cpp/detail/grid.hxx" + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace bioimage_cpp::skeleton::detail { + +struct BallInvalidationEntry { + double distance = 0.0; + std::size_t source = 0; + std::size_t voxel = 0; +}; + +struct BallInvalidationGreater { + bool operator()( + const BallInvalidationEntry &first, + const BallInvalidationEntry &second + ) const noexcept { + if (first.distance != second.distance) { + return first.distance > second.distance; + } + if (first.source != second.source) { + return first.source > second.source; + } + return first.voxel > second.voxel; + } +}; + +// Invalidate the active, 26-connected part of each strict physical-radius ball. +// Inactive voxels are barriers. Active source voxels are always invalidated. +inline std::size_t invalidate_path_balls( + std::vector &active, + const std::span path, + const std::span radii, + const std::vector &shape, + const std::array &spacing +) { + if (shape.size() != 3) { + throw std::invalid_argument("ball invalidation requires a 3D shape"); + } + if (path.size() != radii.size()) { + throw std::invalid_argument("ball invalidation path and radii must match"); + } + if (active.size() != bioimage_cpp::detail::number_of_elements(shape)) { + throw std::invalid_argument("ball invalidation active mask shape mismatch"); + } + + const auto strides = bioimage_cpp::detail::c_order_strides(shape); + std::vector> source_coordinates(path.size()); + std::priority_queue< + BallInvalidationEntry, + std::vector, + BallInvalidationGreater + > queue; + for (std::size_t source = 0; source < path.size(); ++source) { + if (path[source] >= active.size()) { + throw std::invalid_argument("ball invalidation path voxel is out of bounds"); + } + bioimage_cpp::detail::coords_from_index( + static_cast(path[source]), + strides, + 3, + source_coordinates[source].data() + ); + queue.push({0.0, source, path[source]}); + } + + std::size_t invalidated = 0; + std::array coordinate{}; + while (!queue.empty()) { + const auto entry = queue.top(); + queue.pop(); + if (active[entry.voxel] == 0) { + continue; + } + + active[entry.voxel] = 0; + ++invalidated; + bioimage_cpp::detail::coords_from_index( + static_cast(entry.voxel), + strides, + 3, + coordinate.data() + ); + for (std::ptrdiff_t dz = -1; dz <= 1; ++dz) { + for (std::ptrdiff_t dy = -1; dy <= 1; ++dy) { + for (std::ptrdiff_t dx = -1; dx <= 1; ++dx) { + if (dz == 0 && dy == 0 && dx == 0) { + continue; + } + const std::array neighbor_coordinate{ + coordinate[0] + dz, + coordinate[1] + dy, + coordinate[2] + dx, + }; + if ( + neighbor_coordinate[0] < 0 || + neighbor_coordinate[1] < 0 || + neighbor_coordinate[2] < 0 || + neighbor_coordinate[0] >= shape[0] || + neighbor_coordinate[1] >= shape[1] || + neighbor_coordinate[2] >= shape[2] + ) { + continue; + } + const auto neighbor = static_cast( + neighbor_coordinate[0] * strides[0] + + neighbor_coordinate[1] * strides[1] + + neighbor_coordinate[2] + ); + if (active[neighbor] == 0) { + continue; + } + + double distance_squared = 0.0; + for (std::size_t axis = 0; axis < 3; ++axis) { + const auto delta = static_cast( + neighbor_coordinate[axis] - + source_coordinates[entry.source][axis] + ) * spacing[axis]; + distance_squared += delta * delta; + } + const double distance = std::sqrt(distance_squared); + if (distance < radii[entry.source]) { + queue.push({distance, entry.source, neighbor}); + } + } + } + } + } + return invalidated; +} + +} // namespace bioimage_cpp::skeleton::detail diff --git a/include/bioimage_cpp/skeleton/teasar.hxx b/include/bioimage_cpp/skeleton/teasar.hxx index cfc952a..d8e9ae6 100644 --- a/include/bioimage_cpp/skeleton/teasar.hxx +++ b/include/bioimage_cpp/skeleton/teasar.hxx @@ -8,6 +8,7 @@ #include "bioimage_cpp/distance/grid_dijkstra.hxx" #include "bioimage_cpp/skeleton/detail/compact_grid_dijkstra.hxx" #include "bioimage_cpp/skeleton/detail/components.hxx" +#include "bioimage_cpp/skeleton/detail/invalidation.hxx" #include "bioimage_cpp/skeleton/detail/row_interval_union.hxx" #include @@ -30,6 +31,11 @@ namespace bioimage_cpp::skeleton { +enum class TeasarInvalidation { + Cube, + Ball, +}; + struct TeasarOptions { std::array spacing{1.0, 1.0, 1.0}; double scale = 1.5; @@ -37,6 +43,8 @@ struct TeasarOptions { double pdrf_scale = 100000.0; double pdrf_exponent = 4.0; std::size_t number_of_threads = 1; + TeasarInvalidation invalidation = TeasarInvalidation::Cube; + bool fix_branching = true; }; using VoxelCoordinate = std::array; @@ -102,6 +110,12 @@ inline void validate_options( if (!(std::isfinite(options.pdrf_exponent) && options.pdrf_exponent > 0.0)) { throw std::invalid_argument("pdrf_exponent must be positive and finite"); } + if ( + options.invalidation != TeasarInvalidation::Cube && + options.invalidation != TeasarInvalidation::Ball + ) { + throw std::invalid_argument("invalid TEASAR invalidation mode"); + } } inline std::size_t farthest_foreground( @@ -344,23 +358,52 @@ inline LatticeSkeletonGraph teasar_dense_impl( graph.radii.push_back(dbf[voxel]); vertex_of_voxel[voxel] = static_cast(vertex_id); skeleton_voxels.push_back(voxel); - pdrf[voxel] = 0.0; + if (options.fix_branching) { + pdrf[voxel] = 0.0; + } return vertex_id; }; - add_vertex(root); ConstArrayView padded_view{padded_mask.data(), shape, {}}; ConstArrayView pdrf_view{pdrf.data(), shape, {}}; const distance::DijkstraOptions node_options{ 3, {}, distance::DijkstraCostMode::Node, effective_threads }; + std::vector fixed_predecessors; + if (!options.fix_branching) { + BIOIMAGE_PROFILE_SCOPE(profile, "parental_field") + auto field = distance::dijkstra_distance_field( + padded_view, {root}, node_options, &pdrf_view, true + ); + fixed_predecessors = std::move(field.predecessors); + } + + add_vertex(root); std::vector path; const auto trace_target = [&](const std::size_t target) { - { + if (options.fix_branching) { BIOIMAGE_PROFILE_SCOPE(profile, "path_dijkstra") path = distance::dijkstra_path( padded_view, target, skeleton_voxels, node_options, &pdrf_view ); + } else { + BIOIMAGE_PROFILE_SCOPE(profile, "path_from_parents") + path.clear(); + auto voxel = target; + while (vertex_of_voxel[voxel] < 0) { + path.push_back(voxel); + const auto predecessor = fixed_predecessors[voxel]; + if (predecessor < 0) { + throw std::runtime_error( + "invalid fixed-parent predecessor chain" + ); + } + voxel = static_cast(predecessor); + if (path.size() > n) { + throw std::runtime_error("cycle in fixed-parent predecessor chain"); + } + } + path.push_back(voxel); } std::uint64_t previous = static_cast( @@ -380,35 +423,63 @@ inline LatticeSkeletonGraph teasar_dense_impl( { BIOIMAGE_PROFILE_SCOPE(profile, "invalidation") - for (const auto voxel : path) { - pdrf[voxel] = 0.0; - bioimage_cpp::detail::coords_from_index( - static_cast(voxel), strides, 3, coords.data() - ); - const double radius = - options.scale * static_cast(dbf[voxel]) + options.constant; - if (!std::isfinite(radius)) { - throw std::runtime_error("TEASAR invalidation radius overflowed"); + if (options.invalidation == TeasarInvalidation::Ball) { + std::vector radii; + radii.reserve(path.size()); + for (const auto voxel : path) { + const double radius = + options.scale * static_cast(dbf[voxel]) + + options.constant; + if (!std::isfinite(radius)) { + throw std::runtime_error("TEASAR invalidation radius overflowed"); + } + radii.push_back(radius); } - std::array lo{}; - std::array hi{}; - detail_teasar::invalidation_bounds( - coords, radius, options.spacing, shape, lo, hi + const auto invalidated = detail::invalidate_path_balls( + active, path, radii, shape, options.spacing ); - for (std::ptrdiff_t z = lo[0]; z <= hi[0]; ++z) { - for (std::ptrdiff_t y = lo[1]; y <= hi[1]; ++y) { - for (std::ptrdiff_t x = lo[2]; x <= hi[2]; ++x) { - const auto index = static_cast( - z * strides[0] + y * strides[1] + x - ); - if (active[index] != 0) { - active[index] = 0; - --active_count; + if (invalidated > active_count) { + throw std::runtime_error( + "TEASAR ball invalidation count is inconsistent" + ); + } + active_count -= invalidated; + } else { + for (const auto voxel : path) { + bioimage_cpp::detail::coords_from_index( + static_cast(voxel), strides, 3, coords.data() + ); + const double radius = + options.scale * static_cast(dbf[voxel]) + + options.constant; + if (!std::isfinite(radius)) { + throw std::runtime_error("TEASAR invalidation radius overflowed"); + } + std::array lo{}; + std::array hi{}; + detail_teasar::invalidation_bounds( + coords, radius, options.spacing, shape, lo, hi + ); + for (std::ptrdiff_t z = lo[0]; z <= hi[0]; ++z) { + for (std::ptrdiff_t y = lo[1]; y <= hi[1]; ++y) { + for (std::ptrdiff_t x = lo[2]; x <= hi[2]; ++x) { + const auto index = static_cast( + z * strides[0] + y * strides[1] + x + ); + if (active[index] != 0) { + active[index] = 0; + --active_count; + } } } } } } + if (options.fix_branching) { + for (const auto voxel : path) { + pdrf[voxel] = 0.0; + } + } } }; @@ -708,21 +779,43 @@ inline LatticeSkeletonGraph teasar_compact_impl( graph.radii.push_back(compact_dbf[node]); vertex_of_node[node] = static_cast(vertex_id); skeleton_nodes.push_back(node); - pdrf[node] = Distance{0}; + if (options.fix_branching) { + pdrf[node] = Distance{0}; + } return vertex_id; }; + std::vector fixed_predecessors; + if (!options.fix_branching) { + BIOIMAGE_PROFILE_SCOPE(profile, "parental_field") + detail::compact_node_cost_parental_field( + domain, root, pdrf, dijkstra_workspace, fixed_predecessors + ); + } + add_vertex(root); std::vector path; detail::RowIntervalUnion invalidated_rows( n / static_cast(shape[2]), shape[2] ); const auto trace_target = [&](const std::uint32_t target) { - { + if (options.fix_branching) { BIOIMAGE_PROFILE_SCOPE(profile, "path_dijkstra") detail::compact_node_cost_path( domain, target, skeleton_nodes, pdrf, dijkstra_workspace, path ); + } else { + BIOIMAGE_PROFILE_SCOPE(profile, "path_from_parents") + path.clear(); + auto node = target; + while (vertex_of_node[node] < 0) { + path.push_back(node); + node = fixed_predecessors[node]; + if (path.size() > domain.size()) { + throw std::runtime_error("cycle in fixed-parent predecessor chain"); + } + } + path.push_back(node); } std::uint64_t previous = static_cast( @@ -742,49 +835,83 @@ inline LatticeSkeletonGraph teasar_compact_impl( { BIOIMAGE_PROFILE_SCOPE(profile, "invalidation") - for (const auto node : path) { - pdrf[node] = Distance{0}; - const auto voxel = static_cast(domain.compact_to_full[node]); - bioimage_cpp::detail::coords_from_index( - static_cast(voxel), strides, 3, coords.data() - ); - const double radius = - options.scale * static_cast(compact_dbf[node]) + - options.constant; - if (!std::isfinite(radius)) { - throw std::runtime_error("TEASAR invalidation radius overflowed"); + if (options.invalidation == TeasarInvalidation::Ball) { + std::vector full_path; + std::vector radii; + full_path.reserve(path.size()); + radii.reserve(path.size()); + for (const auto node : path) { + full_path.push_back( + static_cast(domain.compact_to_full[node]) + ); + const double radius = + options.scale * static_cast(compact_dbf[node]) + + options.constant; + if (!std::isfinite(radius)) { + throw std::runtime_error("TEASAR invalidation radius overflowed"); + } + radii.push_back(radius); } - std::array lo{}; - std::array hi{}; - detail_teasar::invalidation_bounds( - coords, radius, options.spacing, shape, lo, hi + const auto invalidated = detail::invalidate_path_balls( + active, full_path, radii, shape, options.spacing ); - for (std::ptrdiff_t z = lo[0]; z <= hi[0]; ++z) { - for (std::ptrdiff_t y = lo[1]; y <= hi[1]; ++y) { - const auto row = static_cast( - z * shape[1] + y - ); - const auto row_begin = static_cast( - z * strides[0] + y * strides[1] - ); - invalidated_rows.insert( - row, - lo[2], - hi[2], - [&](const std::ptrdiff_t begin, const std::ptrdiff_t end) { - for (auto x = begin; x <= end; ++x) { - const auto index = row_begin + - static_cast(x); - if (active[index] != 0) { - active[index] = 0; - --active_count; + if (invalidated > active_count) { + throw std::runtime_error( + "TEASAR ball invalidation count is inconsistent" + ); + } + active_count -= invalidated; + } else { + for (const auto node : path) { + const auto voxel = static_cast( + domain.compact_to_full[node] + ); + bioimage_cpp::detail::coords_from_index( + static_cast(voxel), strides, 3, coords.data() + ); + const double radius = + options.scale * static_cast(compact_dbf[node]) + + options.constant; + if (!std::isfinite(radius)) { + throw std::runtime_error("TEASAR invalidation radius overflowed"); + } + std::array lo{}; + std::array hi{}; + detail_teasar::invalidation_bounds( + coords, radius, options.spacing, shape, lo, hi + ); + for (std::ptrdiff_t z = lo[0]; z <= hi[0]; ++z) { + for (std::ptrdiff_t y = lo[1]; y <= hi[1]; ++y) { + const auto row = static_cast( + z * shape[1] + y + ); + const auto row_begin = static_cast( + z * strides[0] + y * strides[1] + ); + invalidated_rows.insert( + row, + lo[2], + hi[2], + [&](const std::ptrdiff_t begin, const std::ptrdiff_t end) { + for (auto x = begin; x <= end; ++x) { + const auto index = row_begin + + static_cast(x); + if (active[index] != 0) { + active[index] = 0; + --active_count; + } } } - } - ); + ); + } } } } + if (options.fix_branching) { + for (const auto node : path) { + pdrf[node] = Distance{0}; + } + } } }; diff --git a/src/bindings/skeleton.cxx b/src/bindings/skeleton.cxx index c630121..7ed71d0 100644 --- a/src/bindings/skeleton.cxx +++ b/src/bindings/skeleton.cxx @@ -53,6 +53,8 @@ nb::tuple teasar_uint8_impl( const double constant, const double pdrf_scale, const double pdrf_exponent, + const bool ball_invalidation, + const bool fix_branching, const skeleton::TeasarBackend backend, const std::size_t n_threads ) { @@ -82,6 +84,9 @@ nb::tuple teasar_uint8_impl( pdrf_scale, pdrf_exponent, n_threads, + ball_invalidation ? skeleton::TeasarInvalidation::Ball + : skeleton::TeasarInvalidation::Cube, + fix_branching, }; if (backend == skeleton::TeasarBackend::Auto) { result = skeleton::teasar(mask_view, options); @@ -99,11 +104,14 @@ nb::tuple teasar_uint8( const double constant, const double pdrf_scale, const double pdrf_exponent, + const bool ball_invalidation, + const bool fix_branching, const std::size_t n_threads ) { return teasar_uint8_impl( mask, spacing, scale, constant, pdrf_scale, pdrf_exponent, - skeleton::TeasarBackend::Auto, n_threads + ball_invalidation, fix_branching, skeleton::TeasarBackend::Auto, + n_threads ); } @@ -114,6 +122,8 @@ nb::tuple teasar_uint8_backend( const double constant, const double pdrf_scale, const double pdrf_exponent, + const bool ball_invalidation, + const bool fix_branching, const std::string &backend, const std::size_t n_threads ) { @@ -128,8 +138,8 @@ nb::tuple teasar_uint8_backend( throw std::invalid_argument("unknown TEASAR development backend: " + backend); } return teasar_uint8_impl( - mask, spacing, scale, constant, pdrf_scale, pdrf_exponent, selected, - n_threads + mask, spacing, scale, constant, pdrf_scale, pdrf_exponent, + ball_invalidation, fix_branching, selected, n_threads ); } @@ -142,6 +152,8 @@ nb::dict teasar_labels_impl( const double constant, const double pdrf_scale, const double pdrf_exponent, + const bool ball_invalidation, + const bool fix_branching, const std::size_t n_threads ) { if (labels.ndim() != 3) { @@ -172,7 +184,10 @@ nb::dict teasar_labels_impl( constant, pdrf_scale, pdrf_exponent, - n_threads} + n_threads, + ball_invalidation ? skeleton::TeasarInvalidation::Ball + : skeleton::TeasarInvalidation::Cube, + fix_branching} ); } @@ -201,6 +216,8 @@ void bind_skeleton(nb::module_ &m) { nb::arg("constant"), nb::arg("pdrf_scale"), nb::arg("pdrf_exponent"), + nb::arg("ball_invalidation"), + nb::arg("fix_branching"), nb::arg("n_threads"), "Core binary 3D TEASAR skeletonization." ); @@ -213,6 +230,8 @@ void bind_skeleton(nb::module_ &m) { nb::arg("constant"), nb::arg("pdrf_scale"), nb::arg("pdrf_exponent"), + nb::arg("ball_invalidation"), + nb::arg("fix_branching"), nb::arg("backend"), nb::arg("n_threads") = 1, "Development-only TEASAR backend selector." @@ -229,6 +248,8 @@ void bind_skeleton(nb::module_ &m) { nb::arg("constant"), \ nb::arg("pdrf_scale"), \ nb::arg("pdrf_exponent"), \ + nb::arg("ball_invalidation"), \ + nb::arg("fix_branching"), \ nb::arg("n_threads"), \ "Core multi-label 3D TEASAR skeletonization." \ ) diff --git a/src/bindings/skeleton_distributed.cxx b/src/bindings/skeleton_distributed.cxx index ca05585..7b4ab60 100644 --- a/src/bindings/skeleton_distributed.cxx +++ b/src/bindings/skeleton_distributed.cxx @@ -201,12 +201,17 @@ skeleton::TeasarOptions teasar_options( const double constant, const double pdrf_scale, const double pdrf_exponent, + const bool ball_invalidation, + const bool fix_branching, const std::size_t number_of_threads ) { const auto values = spacing_array(spacing); return { values, scale, constant, pdrf_scale, pdrf_exponent, - number_of_threads + number_of_threads, + ball_invalidation ? skeleton::TeasarInvalidation::Ball + : skeleton::TeasarInvalidation::Cube, + fix_branching }; } @@ -285,6 +290,8 @@ nb::tuple block_teasar_uint8( const double constant, const double pdrf_scale, const double pdrf_exponent, + const bool ball_invalidation, + const bool fix_branching, const std::size_t number_of_threads ) { const auto shape = array_shape(mask); @@ -301,7 +308,7 @@ nb::tuple block_teasar_uint8( open_face_policy(open_axes, open_high), teasar_options( spacing, scale, constant, pdrf_scale, pdrf_exponent, - number_of_threads + ball_invalidation, fix_branching, number_of_threads ) ); } @@ -321,6 +328,8 @@ nb::dict block_teasar_labels_t( const double constant, const double pdrf_scale, const double pdrf_exponent, + const bool ball_invalidation, + const bool fix_branching, const std::size_t number_of_threads ) { const auto shape = array_shape(labels); @@ -345,7 +354,7 @@ nb::dict block_teasar_labels_t( open_face_policy(open_axes, open_high), teasar_options( spacing, scale, constant, pdrf_scale, pdrf_exponent, - number_of_threads + ball_invalidation, fix_branching, number_of_threads ) ); } @@ -489,7 +498,9 @@ void bind_skeleton_distributed(nb::module_ &m) { nb::arg("mask"), nb::arg("required_targets"), nb::arg("open_axes"), nb::arg("open_high"), nb::arg("origin"), nb::arg("spacing"), nb::arg("scale"), nb::arg("constant"), - nb::arg("pdrf_scale"), nb::arg("pdrf_exponent"), nb::arg("n_threads") + nb::arg("pdrf_scale"), nb::arg("pdrf_exponent"), + nb::arg("ball_invalidation"), nb::arg("fix_branching"), + nb::arg("n_threads") ); #define BIC_BIND_BLOCK_LABELS(name, type) \ @@ -508,6 +519,7 @@ void bind_skeleton_distributed(nb::module_ &m) { nb::arg("open_high"), nb::arg("origin"), \ nb::arg("spacing"), nb::arg("scale"), nb::arg("constant"), \ nb::arg("pdrf_scale"), nb::arg("pdrf_exponent"), \ + nb::arg("ball_invalidation"), nb::arg("fix_branching"), \ nb::arg("n_threads") \ ) diff --git a/src/bioimage_cpp/skeleton/__init__.py b/src/bioimage_cpp/skeleton/__init__.py index 9a3609e..b8e9501 100644 --- a/src/bioimage_cpp/skeleton/__init__.py +++ b/src/bioimage_cpp/skeleton/__init__.py @@ -22,6 +22,11 @@ np.dtype("int64"): _core._teasar_labels_int64, } +_TEASAR_INVALIDATIONS = { + "cube": False, + "ball": True, +} + def _finite_parameter(value, name: str, *, positive: bool) -> float: try: @@ -42,8 +47,10 @@ def _normalize_teasar_options( constant: float, pdrf_scale: float, pdrf_exponent: float, + invalidation: str, + fix_branching: bool, number_of_threads: int, -) -> tuple[list[float], float, float, float, float, int]: +) -> tuple[list[float], float, float, float, float, bool, bool, int]: spacing_values = _normalize_sampling(spacing, 3, function, name="spacing") scale_value = _finite_parameter(scale, "scale", positive=False) constant_value = _finite_parameter(constant, "constant", positive=False) @@ -51,6 +58,14 @@ def _normalize_teasar_options( pdrf_exponent_value = _finite_parameter( pdrf_exponent, "pdrf_exponent", positive=True ) + try: + ball_invalidation = _TEASAR_INVALIDATIONS[invalidation] + except (KeyError, TypeError) as error: + supported = ", ".join(repr(value) for value in _TEASAR_INVALIDATIONS) + raise ValueError( + f"{function}: invalidation must be one of ({supported}), " + f"got {invalidation!r}" + ) from error n_threads = _normalize_threads(number_of_threads, function) return ( spacing_values, @@ -58,6 +73,8 @@ def _normalize_teasar_options( constant_value, pdrf_scale_value, pdrf_exponent_value, + ball_invalidation, + bool(fix_branching), n_threads, ) @@ -70,13 +87,15 @@ def teasar( constant: float = 0.0, pdrf_scale: float = 100000.0, pdrf_exponent: float = 4.0, + invalidation: str = "cube", + fix_branching: bool = True, number_of_threads: int = 1, ) -> tuple[np.ndarray, np.ndarray, np.ndarray]: """Skeletonize a binary volume with 3D TEASAR. This is a correctness-first implementation of the core TEASAR procedure: distance-to-boundary and distance-from-root fields guide repeated penalized - Dijkstra paths, and a rolling physical invalidation cube determines when + Dijkstra paths, and rolling physical invalidation regions determine when the object has been covered. Production kimimaro heuristics such as soma handling, border stitching, hole filling, and manual targets are not part of this function. @@ -95,6 +114,13 @@ def teasar( pdrf_scale, pdrf_exponent: Scale and exponent of the boundary-avoidance term in the penalized distance-from-root field. + invalidation: + ``"cube"`` uses the default physical axis-aligned invalidation box. + ``"ball"`` invalidates the active, foreground-connected part of each + strict physical-radius ball. + fix_branching: + If true, each path is routed to the existing zero-cost skeleton. If + false, all paths follow one parental field computed from the root. number_of_threads: Thread budget for the exact distance transform. Compact Dijkstra root and rail solves remain sequential because their wavefronts benchmarked @@ -125,7 +151,7 @@ def teasar( raise ValueError(f"{function}: mask must have ndim 3, got ndim={binary.ndim}") options = _normalize_teasar_options( function, spacing, scale, constant, pdrf_scale, pdrf_exponent, - number_of_threads + invalidation, fix_branching, number_of_threads ) return _core._teasar_uint8( binary, @@ -142,6 +168,8 @@ def teasar_labels( constant: float = 0.0, pdrf_scale: float = 100000.0, pdrf_exponent: float = 4.0, + invalidation: str = "cube", + fix_branching: bool = True, number_of_threads: int = 1, ) -> dict[int, tuple[np.ndarray, np.ndarray, np.ndarray]]: """Skeletonize every semantic label in a 3D integer volume. @@ -160,7 +188,8 @@ def teasar_labels( background: Integer value excluded from skeletonization. It must fit the input dtype. Defaults to ``0``. - spacing, scale, constant, pdrf_scale, pdrf_exponent, number_of_threads: + spacing, scale, constant, pdrf_scale, pdrf_exponent, invalidation, + fix_branching, number_of_threads: The same TEASAR parameters and shared thread budget as :func:`teasar`. Returns @@ -198,7 +227,7 @@ def teasar_labels( labels_c = np.ascontiguousarray(labels_array) options = _normalize_teasar_options( function, spacing, scale, constant, pdrf_scale, pdrf_exponent, - number_of_threads + invalidation, fix_branching, number_of_threads ) return run(labels_c, labels_array.dtype.type(background_value), *options) diff --git a/src/bioimage_cpp/skeleton/distributed.py b/src/bioimage_cpp/skeleton/distributed.py index 65ef80a..8f7eac1 100644 --- a/src/bioimage_cpp/skeleton/distributed.py +++ b/src/bioimage_cpp/skeleton/distributed.py @@ -203,6 +203,8 @@ def block_teasar( constant: float = 0.0, pdrf_scale: float = 100000.0, pdrf_exponent: float = 4.0, + invalidation: str = "cube", + fix_branching: bool = True, number_of_threads: int = 1, ) -> SkeletonFragment: """Skeletonize a binary processing block into a global lattice fragment. @@ -212,7 +214,8 @@ def block_teasar( confined to real input voxels. ``required_targets`` contains global coordinates and is normally the union of these faces' targets. One target on an open face becomes a deterministic component root and every other - target is forced onto a rail. + target is forced onto a rail. ``invalidation`` and ``fix_branching`` have + the same behavior as in :func:`bioimage_cpp.skeleton.teasar`. """ function = "block_teasar" binary = _as_binary_input(mask, function) @@ -225,7 +228,7 @@ def block_teasar( open_axes, open_high = _normalize_faces(open_faces) options = _normalize_teasar_options( function, spacing, scale, constant, pdrf_scale, pdrf_exponent, - number_of_threads, + invalidation, fix_branching, number_of_threads, ) return _core._block_teasar_uint8( binary, targets, open_axes, open_high, _normalize_origin(origin), *options @@ -244,6 +247,8 @@ def block_teasar_labels( constant: float = 0.0, pdrf_scale: float = 100000.0, pdrf_exponent: float = 4.0, + invalidation: str = "cube", + fix_branching: bool = True, number_of_threads: int = 1, ) -> dict[int, SkeletonFragment]: """Skeletonize a labeled processing block into global lattice fragments. @@ -252,7 +257,9 @@ def block_teasar_labels( ``required_targets`` maps original labels to global coordinate arrays. A coordinate must contain exactly its mapping key in ``labels``. One open-face target roots each affected component. The result has one global lattice - forest per original non-background label. + forest per original non-background label. ``invalidation`` and + ``fix_branching`` have the same behavior as in + :func:`bioimage_cpp.skeleton.teasar_labels`. """ function = "block_teasar_labels" array = _normalize_labels(labels, function) @@ -274,7 +281,7 @@ def block_teasar_labels( ) options = _normalize_teasar_options( function, spacing, scale, constant, pdrf_scale, pdrf_exponent, - number_of_threads, + invalidation, fix_branching, number_of_threads, ) open_axes, open_high = _normalize_faces(open_faces) return _BLOCK_TEASAR_LABELS[array.dtype]( diff --git a/tests/skeleton/test_distributed.py b/tests/skeleton/test_distributed.py index c9baa9c..85b8373 100644 --- a/tests/skeleton/test_distributed.py +++ b/tests/skeleton/test_distributed.py @@ -79,7 +79,13 @@ def test_border_target_edge_tiebreak_uses_physical_low_edge_distance(): np.testing.assert_array_equal(targets, [[1, 1, 1], [2, 4, 1]]) -def test_block_teasar_required_target_and_existing_teasar_equivalence(): +@pytest.mark.parametrize( + "invalidation, fix_branching", + [("cube", True), ("cube", False), ("ball", True), ("ball", False)], +) +def test_block_teasar_required_target_and_existing_teasar_equivalence( + invalidation, fix_branching +): mask = np.zeros((9, 9, 13), dtype=np.uint8) mask[2:7, 2:7, 1:12] = 1 target = np.array([[106, 206, 310]], dtype=np.int64) @@ -92,10 +98,20 @@ def test_block_teasar_required_target_and_existing_teasar_equivalence(): assert block[0].dtype == np.int64 unconstrained = dist.block_teasar( - mask, open_faces=(), origin=(0, 0, 0), spacing=(2, 1, 0.5) + mask, + open_faces=(), + origin=(0, 0, 0), + spacing=(2, 1, 0.5), + invalidation=invalidation, + fix_branching=fix_branching, ) physical = dist.lattice_to_physical(unconstrained, spacing=(2, 1, 0.5)) - ordinary = bic.skeleton.teasar(mask, spacing=(2, 1, 0.5)) + ordinary = bic.skeleton.teasar( + mask, + spacing=(2, 1, 0.5), + invalidation=invalidation, + fix_branching=fix_branching, + ) for actual, expected in zip(physical, ordinary): np.testing.assert_array_equal(actual, expected) @@ -223,6 +239,30 @@ def test_labeled_open_face_uses_label_specific_distance_boundary(): assert np.all(labels[tuple(graph[0].T)] == -5) +@pytest.mark.parametrize( + "invalidation, fix_branching", + [("cube", True), ("cube", False), ("ball", True), ("ball", False)], +) +def test_block_labels_propagate_invalidation_and_branching_options( + invalidation, fix_branching +): + labels = np.zeros((11, 11, 17), dtype=np.uint16) + labels[2:7, 2:7, 2:7] = 11 + labels[6:11, 6:11, 10:15] = 29 + options = { + "open_faces": (), + "scale": 0.0, + "constant": 2.0, + "invalidation": invalidation, + "fix_branching": fix_branching, + } + result = dist.block_teasar_labels(labels, **options) + for label in (11, 29): + expected = dist.block_teasar(labels == label, **options) + for got, wanted in zip(result[label], expected): + np.testing.assert_array_equal(got, wanted) + + @pytest.mark.parametrize( "targets, message", [ diff --git a/tests/skeleton/test_teasar.py b/tests/skeleton/test_teasar.py index 01a1d3f..245c44c 100644 --- a/tests/skeleton/test_teasar.py +++ b/tests/skeleton/test_teasar.py @@ -6,7 +6,13 @@ def _teasar_backend( - mask, backend, *, spacing=(1.5, 1.0, 1.0), number_of_threads=1 + mask, + backend, + *, + spacing=(1.5, 1.0, 1.0), + ball_invalidation=False, + fix_branching=True, + number_of_threads=1, ): return _core._teasar_uint8_backend( np.ascontiguousarray(mask, dtype=np.uint8), @@ -15,6 +21,8 @@ def _teasar_backend( 1.0, 100000.0, 4.0, + ball_invalidation, + fix_branching, backend, number_of_threads, ) @@ -192,25 +200,98 @@ def test_output_is_deterministic(): np.testing.assert_array_equal(got, expected) -def test_threaded_output_is_deterministic(): +def test_default_options_preserve_cube_invalidation_and_branch_fixing(): + mask = np.zeros((9, 9, 9), dtype=np.uint8) + mask[2:7, 2:7, 2:7] = 1 + default = bic.skeleton.teasar(mask, scale=0.0, constant=2.0) + explicit = bic.skeleton.teasar( + mask, + scale=0.0, + constant=2.0, + invalidation="cube", + fix_branching=True, + ) + for got, expected in zip(default, explicit): + np.testing.assert_array_equal(got, expected) + + +def test_ball_invalidation_uses_a_strict_physical_radius(): + mask = np.zeros((9, 9, 9), dtype=np.uint8) + mask[2:7, 2:7, 2:7] = 1 + cube = bic.skeleton.teasar( + mask, scale=0.0, constant=2.0, invalidation="cube" + ) + ball = bic.skeleton.teasar( + mask, scale=0.0, constant=2.0, invalidation="ball" + ) + _assert_valid_tree(mask, *cube) + _assert_valid_tree(mask, *ball) + assert len(cube[0]) == 5 + assert len(ball[0]) == 23 + + +def test_fix_branching_false_uses_the_fixed_parental_field(): + mask = np.zeros((15, 21, 25), dtype=np.uint8) + mask[7, 10, 2:20] = 1 + mask[7, 3:18, 12] = 1 + mask[3:12, 10, 12] = 1 + fixed = bic.skeleton.teasar(mask, constant=1.0, fix_branching=True) + parental = bic.skeleton.teasar(mask, constant=1.0, fix_branching=False) + _assert_valid_tree(mask, *fixed) + _assert_valid_tree(mask, *parental) + fixed_degree = np.bincount(fixed[1].ravel(), minlength=len(fixed[0])) + parental_degree = np.bincount( + parental[1].ravel(), minlength=len(parental[0]) + ) + assert np.count_nonzero(fixed_degree > 2) == 2 + assert np.count_nonzero(parental_degree > 2) == 1 + + +@pytest.mark.parametrize( + "invalidation, fix_branching", + [("cube", True), ("cube", False), ("ball", True), ("ball", False)], +) +def test_threaded_output_is_deterministic(invalidation, fix_branching): mask = np.zeros((17, 21, 25), dtype=bool) mask[8, 10, 2:18] = True mask[8, 4:17, 17] = True - first = bic.skeleton.teasar(mask, number_of_threads=2) - second = bic.skeleton.teasar(mask, number_of_threads=4) + options = { + "invalidation": invalidation, + "fix_branching": fix_branching, + } + first = bic.skeleton.teasar(mask, number_of_threads=2, **options) + second = bic.skeleton.teasar(mask, number_of_threads=4, **options) for got, expected in zip(first, second): np.testing.assert_array_equal(got, expected) @pytest.mark.parametrize("spacing", [(1.0, 1.0, 1.0), (2.5, 1.25, 0.75)]) -def test_compact_fp64_backends_have_exact_dense_parity(spacing): +@pytest.mark.parametrize( + "ball_invalidation, fix_branching", + [(False, True), (False, False), (True, True), (True, False)], +) +def test_compact_fp64_backends_have_exact_dense_parity( + spacing, ball_invalidation, fix_branching +): zz, yy, xx = np.indices((17, 21, 25)) first = ((zz - 8) ** 2 + (yy - 7) ** 2 <= 3**2) & (xx >= 3) & (xx <= 16) second = ((zz - 8) ** 2 + (xx - 16) ** 2 <= 3**2) & (yy >= 7) & (yy <= 17) mask = first | second - dense = _teasar_backend(mask, "dense-fp64", spacing=spacing) + dense = _teasar_backend( + mask, + "dense-fp64", + spacing=spacing, + ball_invalidation=ball_invalidation, + fix_branching=fix_branching, + ) for backend in ("compact-on-the-fly-fp64", "compact-csr-fp64"): - compact = _teasar_backend(mask, backend, spacing=spacing) + compact = _teasar_backend( + mask, + backend, + spacing=spacing, + ball_invalidation=ball_invalidation, + fix_branching=fix_branching, + ) for got, expected in zip(compact, dense): np.testing.assert_array_equal(got, expected) @@ -285,6 +366,8 @@ def test_many_rail_target_ordering_preserves_dense_parity(): 1.0, 100000.0, 4.0, + False, + True, backend, 1, ) @@ -309,7 +392,7 @@ def test_rejects_non_3d_input(shape): def test_direct_binding_validates_ndim_and_spacing_before_dispatch(): - parameters = (1.5, 1.0, 100000.0, 4.0, 1) + parameters = (1.5, 1.0, 100000.0, 4.0, False, True, 1) with pytest.raises(ValueError, match="mask must have ndim 3, got ndim=2"): _core._teasar_uint8( np.ones((3, 4), dtype=np.uint8), [1.0, 1.0, 1.0], *parameters @@ -329,6 +412,7 @@ def test_direct_binding_validates_ndim_and_spacing_before_dispatch(): ({"constant": np.inf}, "constant"), ({"pdrf_scale": np.nan}, "pdrf_scale"), ({"pdrf_exponent": 0.0}, "pdrf_exponent"), + ({"invalidation": "sphere"}, "invalidation"), ({"number_of_threads": -1}, "number_of_threads"), ], ) diff --git a/tests/skeleton/test_teasar_labels.py b/tests/skeleton/test_teasar_labels.py index 1b526b9..cba8859 100644 --- a/tests/skeleton/test_teasar_labels.py +++ b/tests/skeleton/test_teasar_labels.py @@ -158,6 +158,30 @@ def test_each_label_is_exactly_equal_to_binary_dispatch(threads): np.testing.assert_array_equal(got, wanted) +@pytest.mark.parametrize( + "invalidation, fix_branching", + [("cube", True), ("cube", False), ("ball", True), ("ball", False)], +) +def test_each_label_propagates_invalidation_and_branching_options( + invalidation, fix_branching +): + labels = np.zeros((11, 11, 17), dtype=np.uint16) + labels[2:7, 2:7, 2:7] = 11 + labels[6:11, 6:11, 10:15] = 29 + options = { + "scale": 0.0, + "constant": 2.0, + "invalidation": invalidation, + "fix_branching": fix_branching, + "number_of_threads": 2, + } + result = bic.skeleton.teasar_labels(labels, **options) + for label in (11, 29): + expected = bic.skeleton.teasar(labels == label, **options) + for got, wanted in zip(result[label], expected): + np.testing.assert_array_equal(got, wanted) + + @pytest.mark.parametrize( "dtype, labels_values", [ @@ -258,7 +282,7 @@ def test_rejects_wrong_dimensionality(shape): def test_direct_binding_validates_ndim_and_spacing(): - parameters = (1.5, 1.0, 100000.0, 4.0, 1) + parameters = (1.5, 1.0, 100000.0, 4.0, False, True, 1) with pytest.raises(ValueError, match="labels must have ndim 3"): _core._teasar_labels_uint32( np.zeros((3, 4), np.uint32), 0, [1.0, 1.0, 1.0], *parameters From 10acaddadc21c85c1827206f5b1edcc3d9f1803b Mon Sep 17 00:00:00 2001 From: Constantin Pape Date: Tue, 4 Aug 2026 21:26:33 +0200 Subject: [PATCH 2/3] Optimize new teasar features --- development/skeleton/PERFORMANCE_NOTES.md | 171 ++++- development/skeleton/benchmark_teasar.py | 256 +++++-- development/skeleton/benchmark_teasar_mrc.py | 711 ++++++++++++++++-- include/bioimage_cpp/detail/profile.hxx | 36 +- .../distance/distance_transform.hxx | 67 +- .../skeleton/detail/compact_grid_dijkstra.hxx | 31 + .../skeleton/detail/components.hxx | 3 + .../skeleton/detail/invalidation.hxx | 228 ++++++ include/bioimage_cpp/skeleton/teasar.hxx | 533 +++++++++++-- src/bindings/skeleton.cxx | 88 +++ tests/skeleton/test_teasar.py | 111 +++ 11 files changed, 2031 insertions(+), 204 deletions(-) diff --git a/development/skeleton/PERFORMANCE_NOTES.md b/development/skeleton/PERFORMANCE_NOTES.md index 801452d..c02ebfe 100644 --- a/development/skeleton/PERFORMANCE_NOTES.md +++ b/development/skeleton/PERFORMANCE_NOTES.md @@ -609,6 +609,11 @@ semantic keys, paired binary/label parity, and raw samples before writing JSON. ## Real filament MRC target-selection follow-up +The current comparison of cube/ball invalidation and both branching modes is +in `development/skeleton/TEASAR_OPTIONS_COMPARISON.md`. It includes matched +Kimimaro quality metrics and one-worker/eight-worker timings on the centered +200-voxel crop. + The synthetic branching tubes used above require only a small number of rails, so their repeated full-domain target scan was negligible. A real binary mask with many long, closely packed filaments exposed a different scaling regime. @@ -629,8 +634,8 @@ python development/skeleton/benchmark_teasar_mrc.py \ The harness memory-maps the MRC, takes centered nested crops, copies each crop to `uint8` outside the timed region, and reports shape, foreground count, -component count (`V - E`), graph size, raw samples, median, and minimum. It -compares vertices, edges, and radii array-exactly when several worker counts +topology metrics, graph size, raw samples, median, and minimum. It compares +bioimage-cpp vertices, edges, and radii array-exactly when several worker counts are requested. The default is one no-warmup measurement at eight workers for 12.5%, 25%, and 50% crops; the full volume is opt-in because it was initially more than a three-minute call. @@ -791,3 +796,165 @@ On the final `128^3` calls, the last measured phase breakdowns were: Merge and cycle removal remain negligible. More blocks mainly increase face analysis and Python orchestration. The comparison intentionally does not apply Kimimaro's heuristic nearest-component joining, dust removal, or tick pruning. + +## Ball-invalidation optimization pass (2026-08-04) + +This pass optimized the optional ball invalidation mode without changing its +results. It also reduced repeated distance-transform work in ordinary binary +TEASAR. Cube invalidation and `fix_branching=True` remain the public defaults. + +### Implemented changes + +The compact ball invalidator now stores one best pending entry per compact +foreground voxel. Generation stamps avoid clearing this state between paths. +The queue retains the original `(distance, source, node)` order, and stale heap +entries are rejected when removed. Neighbor metadata and source coordinates +are reused across paths. + +Ordinary binary TEASAR can now compute one shared distance transform before +component dispatch. The automatic strategy selects the shared transform only +when all of these conditions hold: + +- the input has at least two components; +- the shared volume is at most 75% of the sum of the component crop volumes; +- the estimated scratch memory, including a 10% allowance, is at most 256 MiB. + +Label and block APIs retain local component transforms. The distance-transform +public API is unchanged. A private development binding can force local or +shared execution for exactness and benchmark checks. + +The profile utility now separates top-level wall time from summed component +work. It also aggregates queue counters across component workers. These +diagnostics compile to no-ops unless `BIOIMAGE_PROFILE` is enabled. + +### Real-mask result + +The headline input is the centered `200 x 200 x 200` crop from +`examples/skeleton/00004_gt_mask.mrc`. It contains 585,875 foreground voxels +in 72 bioimage-cpp components. See `TEASAR_OPTIONS_COMPARISON.md` for the full +topology comparison and reproduction command. + +The table separates the invalidation-queue change from the shared-transform +change. Each final value is the median of seven calls after two warmups. + +| ball setting | workers | baseline | queue only | final | final change | +| --- | ---: | ---: | ---: | ---: | ---: | +| fix | 1 | 1.985 s | 1.121 s | 0.944 s | -52.4% | +| fix | 8 | 0.518 s | 0.306 s | 0.272 s | -47.5% | +| parental | 1 | 1.945 s | 1.096 s | 0.911 s | -53.1% | +| parental | 8 | 0.516 s | 0.311 s | 0.276 s | -46.6% | + +The queue change reduced baseline time by 39.8--43.7%. The shared transform +then reduced queue-only time by a further 11.2--16.8%. The final one-worker +ball modes are 1.94--1.99x faster than matched Kimimaro 5.8.1 modes. The final +eight-worker ball modes are approximately 3.5x faster. + +Cube modes also benefit from the shared transform. Cube/fix fell from 0.975 s +to 0.713 s with one worker and from 0.283 s to 0.219 s with eight workers. +Cube/parental fell from 0.933 s to 0.692 s and from 0.266 s to 0.210 s. + +### Profile result + +A one-worker profile-build call for ball/fix reported 980.9 ms at the top +level. The main component phase totals were: + +| phase | time | +| --- | ---: | +| shared transform setup | 21.8 ms | +| shared distance transform | 151.4 ms | +| compact-domain construction | 27.0 ms | +| root Dijkstra | 283.6 ms | +| PDRF construction | 15.0 ms | +| target selection | 10.6 ms | +| rail-path Dijkstra | 167.4 ms | +| ball invalidation | 273.8 ms | + +The invalidator received 5,976,824 queue offers. It accepted 693,630 offers, +removed 107,755 stale entries, invalidated 585,875 voxels, and reached a peak +heap size of 18,059. Coalescing rejected 88.4% of offers before heap insertion. +The invalidation phase fell from approximately 1.074 s to 273.8 ms. The shared +transform replaced approximately 423 ms of local component transforms with +151 ms of transform work. + +Root Dijkstra, invalidation, and rail-path Dijkstra are now the three largest +one-worker phases. Further work should start with these measured phases after +this fifth optimization step. + +### Memory result + +Fresh worker processes sampled aggregate RSS for each worker and its children. +The table reports the median incremental peak from three processes after +imports and input allocation. + +| implementation and setting | one worker | eight workers | +| --- | ---: | ---: | +| bioimage-cpp ball/fix | 52.7 MiB | 77.1 MiB | +| bioimage-cpp ball/parental | 52.7 MiB | 77.1 MiB | +| Kimimaro ball/fix | 89.8 MiB | 1,575 MiB | +| Kimimaro ball/parental | 89.8 MiB | 1,555 MiB | + +The previous ball/fix implementation used 42.5 MiB with one worker and +129.8 MiB with eight workers. The shared transform adds approximately 10 MiB +to this one-worker process measurement. It saves approximately 53 MiB with +eight workers by avoiding concurrent local transforms. The selected shared +transform estimated 53.9 MB of scratch memory before execution. + +### Synthetic option matrix + +The development benchmark now accepts `cube-fix`, `cube-parental`, `ball-fix`, +and `ball-parental`. This command runs all settings across binary and label +dispatch cases: + +```bash +python development/skeleton/benchmark_teasar.py \ + --large --suite all \ + --settings cube-fix cube-parental ball-fix ball-parental \ + --threads 1 8 --repeats 5 --warmup 1 \ + --json /tmp/teasar_all_settings.json +``` + +The branching-tube medians from the final normal build were: + +| volume | setting | one worker | eight workers | +| --- | --- | ---: | ---: | +| `128^3` | cube/fix | 44.60 ms | 28.52 ms | +| `128^3` | cube/parental | 41.80 ms | 30.06 ms | +| `128^3` | ball/fix | 52.67 ms | 36.87 ms | +| `128^3` | ball/parental | 50.44 ms | 40.26 ms | +| `192^3` | cube/fix | 125.03 ms | 91.96 ms | +| `192^3` | cube/parental | 128.22 ms | 88.85 ms | +| `192^3` | ball/fix | 149.21 ms | 105.67 ms | +| `192^3` | ball/parental | 153.78 ms | 113.57 ms | +| `256^3` | cube/fix | 338.90 ms | 224.21 ms | +| `256^3` | cube/parental | 363.90 ms | 236.88 ms | +| `256^3` | ball/fix | 397.00 ms | 277.87 ms | +| `256^3` | ball/parental | 417.51 ms | 274.86 ms | + +The packed binary and packed semantic-label inputs contain the same foreground +geometry. They produced array-exact forests for every setting. The 27-label +native medians ranged from 54.06 to 71.92 ms with one worker and from 18.70 to +21.03 ms with eight workers. The imbalanced 65-label case did not scale because +one component dominates its small tasks. This result confirms the existing +component-dispatch limitation across all settings. + +The MRC size sweep also measured both ball modes. The 12.5% crop contains +152,662 foreground voxels; ball/fix took 218 ms with one worker and 101 ms with +eight workers. Ball/parental took 213 ms and 96 ms. The 25% crop contains +972,194 foreground voxels; the corresponding medians were 2.015 s and 1.603 s +for ball/fix, and 1.896 s and 1.491 s for ball/parental. The automatic strategy +selected the shared transform for the 25% crop and estimated 76.3 MB of +scratch memory. The more expensive 50% option matrix was not rerun in this +pass. + +### Correctness gates + +Randomized tests compare dense and compact ball invalidation across seeds, +anisotropic spacing, and both branching modes. Forced local and shared distance +transforms must produce array-exact vertices, edges, and radii. Additional +tests cover automatic selection, the scratch-memory fallback, and the +shared-volume fallback. + +The example MRC skeleton statistics remained unchanged through both +optimization stages. Worker-count checks also remained array-exact. Final +verification used the normal build with profile-only instrumentation disabled: +`1421 passed`. diff --git a/development/skeleton/benchmark_teasar.py b/development/skeleton/benchmark_teasar.py index 232f479..08fca31 100644 --- a/development/skeleton/benchmark_teasar.py +++ b/development/skeleton/benchmark_teasar.py @@ -52,6 +52,20 @@ class Workload: components: int +@dataclass(frozen=True) +class Setting: + invalidation: str + fix_branching: bool + + +SETTINGS = { + "cube-fix": Setting("cube", True), + "cube-parental": Setting("cube", False), + "ball-fix": Setting("ball", True), + "ball-parental": Setting("ball", False), +} + + def draw_ball(mask: np.ndarray, center: np.ndarray, radius: int) -> None: center = np.rint(center).astype(int) lo = np.maximum(center - radius, 0) @@ -245,7 +259,10 @@ def kimimaro_parameters(parameters): } -def kimimaro_call(volume, scenario, spacing, parameters, number_of_threads=1): +def kimimaro_call( + volume, scenario, spacing, parameters, fix_branching, + number_of_threads=1, +): import kimimaro kwargs = {} @@ -257,7 +274,7 @@ def kimimaro_call(volume, scenario, spacing, parameters, number_of_threads=1): anisotropy=spacing, dust_threshold=0, progress=False, - fix_branching=True, + fix_branching=fix_branching, fix_borders=False, fill_holes=False, parallel=number_of_threads, @@ -265,7 +282,9 @@ def kimimaro_call(volume, scenario, spacing, parameters, number_of_threads=1): ) -def python_label_loop(volume, spacing, parameters, number_of_threads): +def python_label_loop( + volume, spacing, parameters, setting, number_of_threads +): output = {} for label in np.unique(volume): if label == 0: @@ -274,12 +293,16 @@ def python_label_loop(volume, spacing, parameters, number_of_threads): volume == label, spacing=spacing, number_of_threads=number_of_threads, + invalidation=setting.invalidation, + fix_branching=setting.fix_branching, **parameters, ) return output -def bic_backend_call(mask, spacing, parameters, backend, number_of_threads=1): +def bic_backend_call( + mask, spacing, parameters, setting, backend, number_of_threads=1 +): """Call a development-only C++ backend without changing the public API.""" return _core._teasar_uint8_backend( mask, @@ -288,8 +311,8 @@ def bic_backend_call(mask, spacing, parameters, backend, number_of_threads=1): parameters["constant"], parameters["pdrf_scale"], parameters["pdrf_exponent"], - False, - True, + setting.invalidation == "ball", + setting.fix_branching, backend, number_of_threads, ) @@ -338,11 +361,12 @@ def validate_result(workload, backend, result, counts): def make_backends(workload, args, spacing, parameters): if args.sequential_backends: + setting = SETTINGS[args.settings[0]] return [ ( backend, lambda mask, backend=backend: bic_backend_call( - mask, spacing, parameters, backend + mask, spacing, parameters, setting, backend ), count_bic, ) @@ -351,40 +375,59 @@ def make_backends(workload, args, spacing, parameters): ) ] output = [] - for threads in args.threads: - if workload.scenario == "binary": - function = lambda volume, threads=threads: bic.skeleton.teasar( - volume, spacing=spacing, number_of_threads=threads, **parameters - ) - counter = count_bic - else: - function = lambda volume, threads=threads: bic.skeleton.teasar_labels( - volume, spacing=spacing, number_of_threads=threads, **parameters - ) - counter = count_label_dict - output.append((f"bioimage-cpp/t{threads}", function, counter)) - if workload.scenario == "labels" and args.python_loop: - output.extend( - ( - f"python-label-loop/t{threads}", - lambda volume, threads=threads: python_label_loop( - volume, spacing, parameters, threads - ), - count_label_dict, + for setting_name in args.settings: + setting = SETTINGS[setting_name] + for threads in args.threads: + if workload.scenario == "binary": + function = lambda volume, threads=threads, setting=setting: ( + bic.skeleton.teasar( + volume, spacing=spacing, number_of_threads=threads, + invalidation=setting.invalidation, + fix_branching=setting.fix_branching, + **parameters, + ) + ) + counter = count_bic + else: + function = lambda volume, threads=threads, setting=setting: ( + bic.skeleton.teasar_labels( + volume, spacing=spacing, number_of_threads=threads, + invalidation=setting.invalidation, + fix_branching=setting.fix_branching, + **parameters, + ) + ) + counter = count_label_dict + output.append(( + f"bioimage-cpp/{setting_name}/t{threads}", function, counter + )) + if workload.scenario == "labels" and args.python_loop: + output.extend( + ( + f"python-label-loop/{setting_name}/t{threads}", + lambda volume, threads=threads, setting=setting: python_label_loop( + volume, spacing, parameters, setting, threads + ), + count_label_dict, + ) + for threads in args.threads ) - for threads in args.threads - ) if args.kimimaro: - output.extend( - ( - f"kimimaro/t{threads}", - lambda volume, threads=threads: kimimaro_call( - volume, workload.scenario, spacing, parameters, threads - ), - count_kimimaro_dict, + branching_modes = { + SETTINGS[name].fix_branching for name in args.settings + } + for fix_branching in sorted(branching_modes, reverse=True): + mode = "fix" if fix_branching else "parental" + output.extend( + ( + f"kimimaro/ball-{mode}/t{threads}", + lambda volume, threads=threads, fix=fix_branching: kimimaro_call( + volume, workload.scenario, spacing, parameters, fix, threads + ), + count_kimimaro_dict, + ) + for threads in args.threads ) - for threads in args.threads - ) return output @@ -456,6 +499,7 @@ def memory_worker(args): "scale": 1.5, "constant": 1.0, "pdrf_scale": 100000.0, "pdrf_exponent": 4.0, } + setting = SETTINGS[args.memory_setting] if args.memory_backend == "kimimaro": importlib.import_module("kimimaro") # import before the baseline sample gc.collect() @@ -467,19 +511,25 @@ def memory_worker(args): if selected.scenario == "binary": result = bic.skeleton.teasar( selected.volume, spacing=spacing, - number_of_threads=args.memory_thread, **parameters + number_of_threads=args.memory_thread, + invalidation=setting.invalidation, + fix_branching=setting.fix_branching, + **parameters, ) counts = count_bic(result) else: result = bic.skeleton.teasar_labels( selected.volume, spacing=spacing, - number_of_threads=args.memory_thread, **parameters + number_of_threads=args.memory_thread, + invalidation=setting.invalidation, + fix_branching=setting.fix_branching, + **parameters, ) counts = count_label_dict(result) else: result = kimimaro_call( selected.volume, selected.scenario, spacing, parameters, - args.memory_thread + setting.fix_branching, args.memory_thread ) counts = count_kimimaro_dict(result) payload = { @@ -487,6 +537,7 @@ def memory_worker(args): "scenario": selected.scenario, "case": selected.name, "threads": args.memory_thread, + "setting": args.memory_setting, "elapsed_s": perf_counter() - start, "input_nbytes": selected.volume.nbytes, "vertices": counts[0], @@ -498,7 +549,9 @@ def memory_worker(args): return 0 -def run_memory_probe(script, workload, backend, threads, tier): +def run_memory_probe( + script, workload, backend, setting, threads, tier, repeat +): command = [ sys.executable, script, "--memory-worker", @@ -506,6 +559,7 @@ def run_memory_probe(script, workload, backend, threads, tier): "--memory-case", workload.name, "--memory-scenario", workload.scenario, "--memory-thread", str(threads), + "--memory-setting", setting, "--tier", tier, ] process = subprocess.Popen( @@ -533,6 +587,7 @@ def run_memory_probe(script, workload, backend, threads, tier): ) peak = max(peak, baseline + self_incremental) payload.update({ + "repeat": repeat, "baseline_process_tree_rss_kib": baseline, "peak_process_tree_rss_kib": peak, "incremental_peak_kib": max(0, peak - baseline), @@ -554,10 +609,13 @@ def run_memory_probes(args, selected_workloads): rows = [] for workload in selected: for backend in (("bioimage-cpp", "kimimaro") if args.kimimaro else ("bioimage-cpp",)): - for threads in args.memory_threads: - rows.append(run_memory_probe( - script, workload, backend, threads, args.tier - )) + for setting in args.settings: + for threads in args.memory_threads: + for repeat in range(args.memory_repeats): + rows.append(run_memory_probe( + script, workload, backend, setting, threads, + args.tier, repeat + )) return rows @@ -578,18 +636,26 @@ def parse_args(): parser.add_argument("--kimimaro", action="store_true") parser.add_argument("--python-loop", action=argparse.BooleanOptionalAction, default=True) parser.add_argument("--sequential-backends", action="store_true") + parser.add_argument( + "--settings", nargs="+", choices=tuple(SETTINGS), default=["cube-fix"] + ) parser.add_argument("--json", default="", help="optional JSON result path") parser.add_argument("--threads", type=int, nargs="+", default=[1]) parser.add_argument("--memory", action="store_true") parser.add_argument("--memory-threads", type=int, nargs="+", default=[1, 4]) + parser.add_argument("--memory-repeats", type=int, default=3) parser.add_argument("--memory-worker", action="store_true", help=argparse.SUPPRESS) parser.add_argument("--memory-backend", choices=("bioimage-cpp", "kimimaro"), help=argparse.SUPPRESS) parser.add_argument("--memory-case", help=argparse.SUPPRESS) parser.add_argument("--memory-scenario", choices=("binary", "labels"), help=argparse.SUPPRESS) parser.add_argument("--memory-thread", type=int, help=argparse.SUPPRESS) + parser.add_argument( + "--memory-setting", choices=tuple(SETTINGS), help=argparse.SUPPRESS + ) args = parser.parse_args() if args.tier is None: args.tier = "small" if args.small else "large" if args.large else "default" + args.settings = list(dict.fromkeys(args.settings)) return args @@ -603,6 +669,10 @@ def main() -> int: raise SystemExit("--threads must contain positive worker counts") if args.sequential_backends and args.suite != "binary": raise SystemExit("--sequential-backends requires --suite binary") + if args.sequential_backends and len(args.settings) != 1: + raise SystemExit("--sequential-backends requires exactly one setting") + if args.memory_repeats < 1: + raise SystemExit("--memory-repeats must be positive") if args.kimimaro and importlib.util.find_spec("kimimaro") is None: raise SystemExit("--kimimaro requested, but kimimaro is not installed") @@ -654,13 +724,21 @@ def main() -> int: foreground = int(np.count_nonzero(workload.volume)) runs = number_of_runs(workload.volume, workload.scenario) if not args.sequential_backends: - reference = results[f"bioimage-cpp/t{args.threads[0]}"] - for threads in args.threads[1:]: - candidate = results[f"bioimage-cpp/t{threads}"] - if not exact_bic_result(reference, candidate, workload.scenario): - raise RuntimeError( - f"{workload.name}: worker count {threads} changed output" - ) + for setting_name in args.settings: + reference = results[ + f"bioimage-cpp/{setting_name}/t{args.threads[0]}" + ] + for threads in args.threads[1:]: + candidate = results[ + f"bioimage-cpp/{setting_name}/t{threads}" + ] + if not exact_bic_result( + reference, candidate, workload.scenario + ): + raise RuntimeError( + f"{workload.name}/{setting_name}: worker count " + f"{threads} changed output" + ) for name, _, count_result in backends: result = results[name] counts = count_result(result) @@ -670,6 +748,13 @@ def main() -> int: "scenario": workload.scenario, "case": workload.name, "backend": name, + "setting": next( + ( + setting_name for setting_name in args.settings + if f"/{setting_name}/" in name + ), + None, + ), "number_of_threads": int(name.rsplit("/t", 1)[1]) if "/t" in name else 1, "shape": list(workload.volume.shape), "full_voxels": int(workload.volume.size), @@ -695,31 +780,46 @@ def main() -> int: if args.kimimaro and not args.sequential_backends: by_name = {row["backend"]: row for row in rows if row["case"] == workload.name} - for threads in args.threads: - bio = by_name[f"bioimage-cpp/t{threads}"]["median_s"] - kimi = by_name[f"kimimaro/t{threads}"]["median_s"] - print(f" bioimage-cpp/t{threads} / kimimaro/t{threads}: {bio / kimi:.2f}x") + for setting_name in args.settings: + mode = "fix" if SETTINGS[setting_name].fix_branching else "parental" + for threads in args.threads: + bio = by_name[ + f"bioimage-cpp/{setting_name}/t{threads}" + ]["median_s"] + kimi = by_name[ + f"kimimaro/ball-{mode}/t{threads}" + ]["median_s"] + print( + f" bioimage-cpp/{setting_name}/t{threads} / " + f"kimimaro/ball-{mode}/t{threads}: {bio / kimi:.2f}x" + ) if args.suite == "all" and not args.sequential_backends: - for threads in args.threads: - binary = stored_results[("packed-binary", f"bioimage-cpp/t{threads}")] - labeled = stored_results[("packed-distinct", f"bioimage-cpp/t{threads}")] - flattened_vertices = [] - flattened_radii = [] - flattened_edges = [] - offset = 0 - for skeleton in labeled.values(): - flattened_vertices.append(skeleton[0]) - flattened_radii.append(skeleton[2]) - flattened_edges.append(skeleton[1] + offset) - offset += len(skeleton[0]) - flattened = ( - np.concatenate(flattened_vertices), - np.concatenate(flattened_edges), - np.concatenate(flattened_radii), - ) - if not all(np.array_equal(a, b) for a, b in zip(binary, flattened)): - raise RuntimeError("paired binary and labeled dispatch lost exact parity") + for setting_name in args.settings: + for threads in args.threads: + backend = f"bioimage-cpp/{setting_name}/t{threads}" + binary = stored_results[("packed-binary", backend)] + labeled = stored_results[("packed-distinct", backend)] + flattened_vertices = [] + flattened_radii = [] + flattened_edges = [] + offset = 0 + for skeleton in labeled.values(): + flattened_vertices.append(skeleton[0]) + flattened_radii.append(skeleton[2]) + flattened_edges.append(skeleton[1] + offset) + offset += len(skeleton[0]) + flattened = ( + np.concatenate(flattened_vertices), + np.concatenate(flattened_edges), + np.concatenate(flattened_radii), + ) + if not all( + np.array_equal(a, b) for a, b in zip(binary, flattened) + ): + raise RuntimeError( + "paired binary and labeled dispatch lost exact parity" + ) print("paired packed binary/multi-label exact parity: True") payload = { @@ -731,6 +831,7 @@ def main() -> int: "warmup": args.warmup, "spacing": spacing, "parameters": parameters, + "settings": args.settings, "results": rows, "memory": memory_rows, } @@ -739,7 +840,8 @@ def main() -> int: for row in memory_rows: print( f" {row['scenario']:>6} {row['case']:>22} " - f"{row['backend']:>12}/t{row['threads']}: " + f"{row['backend']:>12}/{row['setting']}/t{row['threads']} " + f"run {row['repeat'] + 1}: " f"{row['incremental_peak_kib'] / 1024:.1f} MiB incremental" ) if args.json: diff --git a/development/skeleton/benchmark_teasar_mrc.py b/development/skeleton/benchmark_teasar_mrc.py index ef1e8cc..30af83c 100644 --- a/development/skeleton/benchmark_teasar_mrc.py +++ b/development/skeleton/benchmark_teasar_mrc.py @@ -1,23 +1,98 @@ +"""Benchmark TEASAR settings on a real binary MRC volume. + +The benchmark supports the four bioimage-cpp combinations of invalidation +geometry and branching behavior. It can also run both Kimimaro branching +modes. Timing excludes MRC loading, crop extraction, and graph statistics. + +Examples +-------- +Run the comparison used by ``TEASAR_OPTIONS_COMPARISON.md``:: + + python development/skeleton/benchmark_teasar_mrc.py \ + --crop-size 200 --pixel-size 10 --scale 0 --constant 70 \ + --ball-constant 140 --all-settings --kimimaro \ + --threads 1 8 --repeats 7 --warmup 2 \ + --json /tmp/teasar_options_mrc.json + +Run the original size sweep with the default TEASAR behavior:: + + python development/skeleton/benchmark_teasar_mrc.py +""" + from __future__ import annotations import argparse +from collections.abc import Callable +from dataclasses import dataclass +from functools import partial +import gc import importlib.metadata +import importlib.util import json import os from pathlib import Path import platform import random from statistics import median +import subprocess import sys +import time from time import perf_counter import mrcfile import numpy as np import bioimage_cpp as bic +from bioimage_cpp import _core +from benchmark_teasar import proc_tree_rss_kib + +try: + import resource +except ImportError: # pragma: no cover - unavailable on Windows + resource = None + + +@dataclass(frozen=True) +class Setting: + name: str + implementation: str + invalidation: str + fix_branching: bool + + +@dataclass(frozen=True) +class Backend: + key: str + setting: Setting + number_of_threads: int + constant: float + function: Callable[[np.ndarray], object] + + +BIOIMAGE_SETTINGS = { + "cube-fix": Setting("bioimage-cpp/cube/fix", "bioimage-cpp", "cube", True), + "cube-parental": Setting( + "bioimage-cpp/cube/parental", "bioimage-cpp", "cube", False + ), + "ball-fix": Setting("bioimage-cpp/ball/fix", "bioimage-cpp", "ball", True), + "ball-parental": Setting( + "bioimage-cpp/ball/parental", "bioimage-cpp", "ball", False + ), +} +KIMIMARO_SETTINGS = ( + Setting("kimimaro/ball/fix", "kimimaro", "ball", True), + Setting("kimimaro/ball/parental", "kimimaro", "ball", False), +) -def centered_crop(shape: tuple[int, ...], fraction: float): +MEMORY_SETTINGS = { + **BIOIMAGE_SETTINGS, + "kimimaro-fix": KIMIMARO_SETTINGS[0], + "kimimaro-parental": KIMIMARO_SETTINGS[1], +} + + +def centered_fraction_crop(shape: tuple[int, ...], fraction: float): shape_array = np.asarray(shape, dtype=np.int64) crop_shape = np.maximum(1, np.floor(shape_array * fraction).astype(np.int64)) begin = (shape_array - crop_shape) // 2 @@ -26,6 +101,14 @@ def centered_crop(shape: tuple[int, ...], fraction: float): return slices, tuple(int(value) for value in begin) +def centered_size_crop(shape: tuple[int, ...], size: int): + if any(size > extent for extent in shape): + raise ValueError(f"crop size {size} exceeds source shape {shape}") + begin = tuple((extent - size) // 2 for extent in shape) + slices = tuple(slice(lo, lo + size) for lo in begin) + return slices, begin + + def package_version(name: str): try: return importlib.metadata.version(name) @@ -33,143 +116,639 @@ def package_version(name: str): return None +def environment(): + return { + "python": sys.version, + "platform": platform.platform(), + "cpu_count": os.cpu_count(), + "numpy": np.__version__, + "bioimage_cpp": package_version("bioimage-cpp"), + "kimimaro": package_version("kimimaro"), + "edt": package_version("edt"), + "mrcfile": package_version("mrcfile"), + "thread_environment": { + key: os.environ.get(key) + for key in ( + "OMP_NUM_THREADS", + "OPENBLAS_NUM_THREADS", + "MKL_NUM_THREADS", + "NUMEXPR_NUM_THREADS", + ) + }, + } + + def exact_result(first, second): return all(np.array_equal(a, b) for a, b in zip(first, second)) -def time_threads(mask, threads, scale, repeats, warmup): +def kimimaro_parameters(parameters): + return { + "scale": parameters["scale"], + "const": parameters["constant"], + "pdrf_scale": parameters["pdrf_scale"], + "pdrf_exponent": int(parameters["pdrf_exponent"]), + "soma_detection_threshold": float("inf"), + "soma_acceptance_threshold": float("inf"), + "soma_invalidation_scale": 1.0, + "soma_invalidation_const": 0.0, + } + + +def kimimaro_call( + mask, + spacing, + parameters, + fix_branching, + number_of_threads, +): + import kimimaro + + return kimimaro.skeletonize( + mask, + teasar_params=kimimaro_parameters(parameters), + anisotropy=spacing, + object_ids=[1], + dust_threshold=0, + progress=False, + fix_branching=fix_branching, + fix_borders=False, + fill_holes=False, + parallel=number_of_threads, + ) + + +def bioimage_cpp_call( + mask, + spacing, + parameters, + invalidation, + fix_branching, + number_of_threads, +): + return bic.skeleton.teasar( + mask, + spacing=spacing, + invalidation=invalidation, + fix_branching=fix_branching, + number_of_threads=number_of_threads, + **parameters, + ) + + +def setting_constant(setting: Setting, constant: float, ball_constant: float): + return ball_constant if setting.invalidation == "ball" else constant + + +def make_backends(args, spacing): + selected = [BIOIMAGE_SETTINGS[name] for name in args.settings] + if args.kimimaro: + selected.extend(KIMIMARO_SETTINGS) + + backends = [] + for setting in selected: + constant = setting_constant(setting, args.constant, args.ball_constant) + parameters = { + "scale": args.scale, + "constant": constant, + "pdrf_scale": args.pdrf_scale, + "pdrf_exponent": args.pdrf_exponent, + } + for threads in args.threads: + key = f"{setting.name}/t{threads}" + if setting.implementation == "bioimage-cpp": + function = partial( + bioimage_cpp_call, + spacing=spacing, + parameters=parameters, + invalidation=setting.invalidation, + fix_branching=setting.fix_branching, + number_of_threads=threads, + ) + else: + function = partial( + kimimaro_call, + spacing=spacing, + parameters=parameters, + fix_branching=setting.fix_branching, + number_of_threads=threads, + ) + backends.append(Backend(key, setting, threads, constant, function)) + return backends + + +def time_backends(mask, backends, repeats, warmup): results = {} - samples = {thread: [] for thread in threads} + samples = {backend.key: [] for backend in backends} for _ in range(warmup): - for thread in threads: - results[thread] = bic.skeleton.teasar( - mask, - scale=scale, - number_of_threads=thread, - ) - rng = random.Random(20260715) + for backend in backends: + results[backend.key] = backend.function(mask) + rng = random.Random(20260804) for _ in range(repeats): - order = list(threads) + order = list(backends) rng.shuffle(order) - for thread in order: + for backend in order: start = perf_counter() - results[thread] = bic.skeleton.teasar( - mask, - scale=scale, - number_of_threads=thread, + results[backend.key] = backend.function(mask) + samples[backend.key].append(perf_counter() - start) + return samples, results + + +def flatten_result(result, implementation): + if implementation == "bioimage-cpp": + vertices, edges, _ = result + return np.asarray(vertices, dtype=np.float64), np.asarray(edges, dtype=np.int64) + + skeletons = result.values() if isinstance(result, dict) else (result,) + vertices = [] + edges = [] + offset = 0 + for skeleton in skeletons: + skeleton_vertices = np.asarray(skeleton.vertices, dtype=np.float64) + skeleton_edges = np.asarray(skeleton.edges, dtype=np.int64).reshape(-1, 2) + vertices.append(skeleton_vertices) + edges.append(skeleton_edges + offset) + offset += len(skeleton_vertices) + if not vertices: + return np.empty((0, 3), np.float64), np.empty((0, 2), np.int64) + return np.concatenate(vertices), np.concatenate(edges) + + +def build_adjacency(number_of_vertices, edges): + if not len(edges): + return ( + np.zeros(number_of_vertices + 1, dtype=np.int64), + np.empty(0, dtype=np.int64), + np.zeros(number_of_vertices, dtype=np.int64), + ) + sources = np.concatenate([edges[:, 0], edges[:, 1]]) + targets = np.concatenate([edges[:, 1], edges[:, 0]]) + order = np.argsort(sources, kind="stable") + targets = targets[order] + degrees = np.bincount(sources, minlength=number_of_vertices) + offsets = np.zeros(number_of_vertices + 1, dtype=np.int64) + np.cumsum(degrees, out=offsets[1:]) + return offsets, targets, degrees + + +def component_count(number_of_vertices, edges): + parents = np.arange(number_of_vertices, dtype=np.int64) + + def find(node): + while parents[node] != node: + parents[node] = parents[parents[node]] + node = int(parents[node]) + return node + + for first, second in edges: + first_root = find(int(first)) + second_root = find(int(second)) + if first_root != second_root: + parents[second_root] = first_root + return len({find(node) for node in range(number_of_vertices)}) + + +def walk_arm(node, first, offsets, targets, degrees, vertices): + previous = node + current = int(first) + length = float(np.linalg.norm(vertices[current] - vertices[node])) + while degrees[current] == 2: + neighbors = targets[offsets[current]:offsets[current + 1]] + next_node = int(neighbors[0]) + if next_node == previous: + next_node = int(neighbors[1]) + length += float(np.linalg.norm(vertices[next_node] - vertices[current])) + previous, current = current, next_node + return current, length + + +def graph_statistics(vertices, edges, spur_length): + offsets, targets, degrees = build_adjacency(len(vertices), edges) + if len(edges): + segment_lengths = np.linalg.norm( + vertices[edges[:, 0]] - vertices[edges[:, 1]], axis=1 + ) + else: + segment_lengths = np.empty(0, dtype=np.float64) + degree_three = np.flatnonzero(degrees == 3) + shortest_arms = [] + spurs = 0 + real_junctions = 0 + for node in degree_three: + arms = [ + walk_arm( + int(node), targets[index], offsets, targets, degrees, vertices ) - samples[thread].append(perf_counter() - start) - reference_thread = threads[0] - for thread in threads[1:]: - if not exact_result(results[reference_thread], results[thread]): - raise RuntimeError( - f"thread count {thread} changed the skeleton compared with " - f"thread count {reference_thread}" + for index in range(offsets[node], offsets[node + 1]) + ] + lengths = [length for _, length in arms] + shortest_arms.append(min(lengths)) + if min(lengths) > spur_length: + real_junctions += 1 + if any( + degrees[end] == 1 and length <= spur_length + for end, length in arms + ): + spurs += 1 + + shortest_arms = np.asarray(shortest_arms, dtype=np.float64) + output = { + "vertices": len(vertices), + "edges": len(edges), + "length_physical": float(segment_lengths.sum()), + "components": component_count(len(vertices), edges), + "degree_1": int(np.count_nonzero(degrees == 1)), + "degree_3": len(degree_three), + "degree_4": int(np.count_nonzero(degrees == 4)), + "spurs": spurs, + "real_junctions": real_junctions, + } + for percentile in (25, 50, 75, 90): + output[f"shortest_arm_p{percentile}"] = ( + float(np.percentile(shortest_arms, percentile)) + if shortest_arms.size + else None + ) + output["spur_percent"] = ( + 100.0 * spurs / len(degree_three) if len(degree_three) else None + ) + output["real_junction_percent"] = ( + 100.0 * real_junctions / len(degree_three) + if len(degree_three) + else None + ) + return output + + +def crop_specs(args, source_shape): + specs = [] + if args.crop_size is not None: + slices, origin = centered_size_crop(source_shape, args.crop_size) + specs.append((f"crop-{args.crop_size}", None, slices, origin)) + if args.include_full: + slices, origin = centered_fraction_crop(source_shape, 1.0) + specs.append(("fraction-1.000", 1.0, slices, origin)) + else: + fractions = set(args.fractions) + if args.include_full: + fractions.add(1.0) + for fraction in sorted(fractions): + slices, origin = centered_fraction_crop(source_shape, fraction) + specs.append((f"fraction-{fraction:.3f}", fraction, slices, origin)) + return specs + + +def memory_worker(args): + if resource is None: + raise RuntimeError("memory worker requires the Unix resource module") + setting = MEMORY_SETTINGS[args.memory_setting] + with mrcfile.mmap(args.input, mode="r", permissive=True) as mrc: + source_shape = tuple(int(value) for value in mrc.data.shape) + specs = crop_specs(args, source_shape) + _, _, slices, _ = specs[args.memory_crop_index] + mask = np.array(mrc.data[slices] > 0, dtype=np.uint8, copy=True) + + spacing = (args.pixel_size,) * 3 + constant = setting_constant(setting, args.constant, args.ball_constant) + parameters = { + "scale": args.scale, + "constant": constant, + "pdrf_scale": args.pdrf_scale, + "pdrf_exponent": args.pdrf_exponent, + } + if setting.implementation == "kimimaro": + importlib.import_module("kimimaro") + gc.collect() + self_peak_before = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss + print("READY", flush=True) + sys.stdin.readline() + start = perf_counter() + if setting.implementation == "bioimage-cpp": + if args.memory_edt_strategy == "auto": + result = bioimage_cpp_call( + mask, spacing, parameters, setting.invalidation, + setting.fix_branching, args.memory_thread, ) - return samples, results + else: + result = _core._teasar_uint8_edt_backend( + mask, + spacing, + parameters["scale"], + parameters["constant"], + parameters["pdrf_scale"], + parameters["pdrf_exponent"], + setting.invalidation == "ball", + setting.fix_branching, + args.memory_edt_strategy, + 256 * 1024 * 1024, + args.memory_thread, + )[:3] + vertices = len(result[0]) + edges = len(result[1]) + else: + result = kimimaro_call( + mask, spacing, parameters, setting.fix_branching, + args.memory_thread, + ) + skeletons = result.values() if isinstance(result, dict) else (result,) + vertices = sum(len(skeleton.vertices) for skeleton in skeletons) + edges = sum(len(skeleton.edges) for skeleton in skeletons) + print(json.dumps({ + "setting": setting.name, + "edt_strategy": args.memory_edt_strategy, + "implementation": setting.implementation, + "threads": args.memory_thread, + "elapsed_s": perf_counter() - start, + "input_nbytes": mask.nbytes, + "vertices": vertices, + "edges": edges, + "self_peak_before_kib": self_peak_before, + "self_peak_after_kib": resource.getrusage(resource.RUSAGE_SELF).ru_maxrss, + }), flush=True) + return 0 + + +def run_memory_probe(args, crop_index, crop_name, setting_key, threads, repeat): + command = [ + sys.executable, + os.path.abspath(__file__), + "--memory-worker", + "--memory-setting", setting_key, + "--memory-crop-index", str(crop_index), + "--memory-thread", str(threads), + "--memory-edt-strategy", args.memory_edt_strategy, + "--input", str(args.input), + "--pixel-size", str(args.pixel_size), + "--scale", str(args.scale), + "--constant", str(args.constant), + "--ball-constant", str(args.ball_constant), + "--pdrf-scale", str(args.pdrf_scale), + "--pdrf-exponent", str(args.pdrf_exponent), + ] + if args.crop_size is not None: + command.extend(("--crop-size", str(args.crop_size))) + else: + command.append("--fractions") + command.extend(str(value) for value in args.fractions) + if args.include_full: + command.append("--include-full") + + process = subprocess.Popen( + command, stdin=subprocess.PIPE, stdout=subprocess.PIPE, + stderr=subprocess.PIPE, text=True, + ) + ready = process.stdout.readline().strip() + if ready != "READY": + _, stderr = process.communicate() + raise RuntimeError(f"memory worker failed before READY: {stderr}") + baseline = proc_tree_rss_kib(process.pid) + peak = baseline + process.stdin.write("go\n") + process.stdin.flush() + while process.poll() is None: + peak = max(peak, proc_tree_rss_kib(process.pid)) + time.sleep(0.002) + peak = max(peak, proc_tree_rss_kib(process.pid)) + stdout, stderr = process.communicate() + if process.returncode != 0: + raise RuntimeError(f"memory worker failed: {stderr}") + payload = json.loads(stdout.strip().splitlines()[-1]) + self_incremental = max( + 0, payload["self_peak_after_kib"] - payload["self_peak_before_kib"] + ) + peak = max(peak, baseline + self_incremental) + payload.update({ + "crop": crop_name, + "repeat": repeat, + "baseline_process_tree_rss_kib": baseline, + "peak_process_tree_rss_kib": peak, + "incremental_peak_kib": max(0, peak - baseline), + }) + return payload + + +def run_memory_probes(args, specs): + if platform.system() != "Linux" or not os.path.isdir("/proc"): + raise RuntimeError("--memory requires Linux /proc") + setting_keys = list(args.settings) + if args.kimimaro: + setting_keys.extend(("kimimaro-fix", "kimimaro-parental")) + rows = [] + for crop_index, (crop_name, _, _, _) in enumerate(specs): + for setting_key in setting_keys: + for threads in args.memory_threads: + for repeat in range(args.memory_repeats): + rows.append(run_memory_probe( + args, crop_index, crop_name, setting_key, threads, repeat + )) + return rows def parse_args(): root = Path(__file__).resolve().parents[2] - parser = argparse.ArgumentParser() + parser = argparse.ArgumentParser( + description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter, + ) parser.add_argument( "--input", type=Path, default=root / "examples" / "skeleton" / "00004_gt_mask.mrc", ) - parser.add_argument( - "--fractions", - type=float, - nargs="+", - default=[0.125, 0.25, 0.5], - ) + crops = parser.add_mutually_exclusive_group() + crops.add_argument("--fractions", type=float, nargs="+") + crops.add_argument("--crop-size", type=int) parser.add_argument("--include-full", action="store_true") parser.add_argument("--threads", type=int, nargs="+", default=[8]) parser.add_argument("--repeats", type=int, default=1) parser.add_argument("--warmup", type=int, default=0) + parser.add_argument("--pixel-size", type=float, default=1.0) parser.add_argument("--scale", type=float, default=3.0) + parser.add_argument("--constant", type=float, default=0.0) + parser.add_argument( + "--ball-constant", + type=float, + help="constant for ball modes; defaults to --constant", + ) + parser.add_argument("--pdrf-scale", type=float, default=100000.0) + parser.add_argument("--pdrf-exponent", type=float, default=4.0) + parser.add_argument("--spur-length", type=float, default=200.0) + parser.add_argument( + "--settings", + nargs="+", + choices=tuple(BIOIMAGE_SETTINGS), + help="bioimage-cpp settings; defaults to cube-fix", + ) + parser.add_argument("--all-settings", action="store_true") + parser.add_argument("--kimimaro", action="store_true") + parser.add_argument("--memory", action="store_true") + parser.add_argument("--memory-threads", type=int, nargs="+", default=[1, 8]) + parser.add_argument("--memory-repeats", type=int, default=3) + parser.add_argument( + "--memory-edt-strategy", choices=("auto", "local", "shared"), + default="auto", + ) + parser.add_argument("--memory-worker", action="store_true", help=argparse.SUPPRESS) + parser.add_argument( + "--memory-setting", choices=tuple(MEMORY_SETTINGS), help=argparse.SUPPRESS + ) + parser.add_argument("--memory-crop-index", type=int, help=argparse.SUPPRESS) + parser.add_argument("--memory-thread", type=int, help=argparse.SUPPRESS) parser.add_argument("--json", type=Path) args = parser.parse_args() + if args.repeats < 1 or args.warmup < 0: parser.error("--repeats must be >= 1 and --warmup must be >= 0") if not args.threads or any(thread < 1 for thread in args.threads): parser.error("--threads must contain positive values") - if any(not 0.0 < fraction <= 1.0 for fraction in args.fractions): + if not args.memory_threads or any(thread < 1 for thread in args.memory_threads): + parser.error("--memory-threads must contain positive values") + if args.memory_repeats < 1: + parser.error("--memory-repeats must be positive") + if args.fractions is not None and any( + not 0.0 < fraction <= 1.0 for fraction in args.fractions + ): parser.error("--fractions values must be in (0, 1]") - fractions = set(args.fractions) - if args.include_full: - fractions.add(1.0) - args.fractions = sorted(fractions) + if args.crop_size is not None and args.crop_size < 1: + parser.error("--crop-size must be positive") + if args.all_settings and args.settings is not None: + parser.error("--all-settings and --settings cannot be combined") + if args.kimimaro and importlib.util.find_spec("kimimaro") is None: + parser.error("--kimimaro requested, but kimimaro is not installed") + for name in ( + "pixel_size", + "pdrf_exponent", + "spur_length", + ): + if not np.isfinite(getattr(args, name)) or getattr(args, name) <= 0.0: + parser.error(f"--{name.replace('_', '-')} must be positive and finite") + for name in ("scale", "constant", "pdrf_scale"): + if not np.isfinite(getattr(args, name)) or getattr(args, name) < 0.0: + parser.error(f"--{name.replace('_', '-')} must be finite and non-negative") + if args.ball_constant is None: + args.ball_constant = args.constant + elif not np.isfinite(args.ball_constant) or args.ball_constant < 0.0: + parser.error("--ball-constant must be finite and non-negative") + + if args.fractions is None and args.crop_size is None: + args.fractions = [0.125, 0.25, 0.5] + args.settings = ( + list(BIOIMAGE_SETTINGS) + if args.all_settings + else list(dict.fromkeys(args.settings or ["cube-fix"])) + ) args.threads = list(dict.fromkeys(args.threads)) return args def main() -> int: args = parse_args() + if args.memory_worker: + return memory_worker(args) rows = [] + spacing = (args.pixel_size,) * 3 + backends = make_backends(args, spacing) header = ( - f"{'fraction':>8} {'threads':>7} {'shape':>18} {'foreground':>11} " - f"{'components':>10} {'vertices':>10} {'median s':>10} {'min s':>10}" + f"{'crop':>14} {'backend':>30} {'t':>2} {'const':>7} " + f"{'shape':>17} {'foreground':>10} {'comp':>5} {'vertices':>8} " + f"{'deg3':>5} {'spurs':>5} {'real':>5} {'median s':>9} {'min s':>9}" ) print(header) print("-" * len(header)) with mrcfile.mmap(args.input, mode="r", permissive=True) as mrc: source = mrc.data source_shape = tuple(int(value) for value in source.shape) - for fraction in args.fractions: - slices, origin = centered_crop(source_shape, fraction) - mask = np.array(source[slices], dtype=np.uint8, copy=True) + specs = crop_specs(args, source_shape) + for crop_name, fraction, slices, origin in specs: + mask = np.array(source[slices] > 0, dtype=np.uint8, copy=True) foreground = int(np.count_nonzero(mask)) - samples, results = time_threads( + samples, results = time_backends( mask, - args.threads, - args.scale, + backends, args.repeats, args.warmup, ) - for thread in args.threads: - vertices, edges, _ = results[thread] - components = len(vertices) - len(edges) + + for setting_name in args.settings: + setting = BIOIMAGE_SETTINGS[setting_name] + reference_key = f"{setting.name}/t{args.threads[0]}" + for threads in args.threads[1:]: + candidate_key = f"{setting.name}/t{threads}" + if not exact_result(results[reference_key], results[candidate_key]): + raise RuntimeError( + f"{setting.name}: thread count {threads} changed the skeleton" + ) + + for backend in backends: + vertices, edges = flatten_result( + results[backend.key], backend.setting.implementation + ) + statistics = graph_statistics(vertices, edges, args.spur_length) + samples_s = samples[backend.key] row = { + "crop": crop_name, "fraction": fraction, "origin": list(origin), "shape": list(mask.shape), "full_voxels": int(mask.size), "foreground_voxels": foreground, - "number_of_threads": thread, - "components": components, - "vertices": len(vertices), - "edges": len(edges), - "samples_s": samples[thread], - "median_s": median(samples[thread]), - "min_s": min(samples[thread]), + "backend": backend.setting.name, + "implementation": backend.setting.implementation, + "invalidation": backend.setting.invalidation, + "fix_branching": backend.setting.fix_branching, + "number_of_threads": backend.number_of_threads, + "constant": backend.constant, + "samples_s": samples_s, + "median_s": median(samples_s), + "min_s": min(samples_s), + **statistics, } rows.append(row) print( - f"{fraction:8.3f} {thread:7d} {str(mask.shape):>18} " - f"{foreground:11d} {components:10d} {len(vertices):10d} " - f"{row['median_s']:10.3f} {row['min_s']:10.3f}" + f"{crop_name:>14} {backend.setting.name:>30} " + f"{backend.number_of_threads:2d} {backend.constant:7.1f} " + f"{str(mask.shape):>17} {foreground:10d} " + f"{statistics['components']:5d} {statistics['vertices']:8d} " + f"{statistics['degree_3']:5d} {statistics['spurs']:5d} " + f"{statistics['real_junctions']:5d} " + f"{row['median_s']:9.3f} {row['min_s']:9.3f}" ) del results del mask + + memory_rows = run_memory_probes(args, specs) if args.memory else [] + if memory_rows: + print("\npeak process-tree RSS") + for row in memory_rows: + print( + f" {row['crop']:>14} {row['setting']:>30}/t{row['threads']} " + f"run {row['repeat'] + 1}: " + f"{row['incremental_peak_kib'] / 1024:.1f} MiB incremental" + ) + if args.json is not None: payload = { - "environment": { - "python": sys.version, - "platform": platform.platform(), - "cpu_count": os.cpu_count(), - "numpy": np.__version__, - "bioimage_cpp": package_version("bioimage-cpp"), - "mrcfile": package_version("mrcfile"), - }, + "environment": environment(), "input": str(args.input.resolve()), "source_shape": list(source_shape), + "crop_size": args.crop_size, "fractions": args.fractions, "threads": args.threads, "repeats": args.repeats, "warmup": args.warmup, - "scale": args.scale, + "spacing": list(spacing), + "parameters": { + "scale": args.scale, + "cube_constant": args.constant, + "ball_constant": args.ball_constant, + "pdrf_scale": args.pdrf_scale, + "pdrf_exponent": args.pdrf_exponent, + "spur_length": args.spur_length, + }, + "settings": args.settings, + "kimimaro": args.kimimaro, "results": rows, + "memory": memory_rows, } args.json.write_text(json.dumps(payload, indent=2), encoding="utf-8") print(f"wrote {args.json}", file=sys.stderr) diff --git a/include/bioimage_cpp/detail/profile.hxx b/include/bioimage_cpp/detail/profile.hxx index 6f1f7ac..1654837 100644 --- a/include/bioimage_cpp/detail/profile.hxx +++ b/include/bioimage_cpp/detail/profile.hxx @@ -8,7 +8,10 @@ // `NullProfiler` is always available; helper templates can request "no // profiling here" via the same type regardless of build mode. namespace bioimage_cpp::detail { -struct NullProfiler {}; +struct NullProfiler { + void merge(const NullProfiler &) noexcept {} + void report(const char * = "[bioimage profile]") const noexcept {} +}; } // namespace bioimage_cpp::detail #ifdef BIOIMAGE_PROFILE @@ -35,8 +38,21 @@ public: } } - void report() const { - std::fprintf(stderr, "[bioimage profile]\n"); + void merge(const Profiler &other) { + for (const auto *name : other.order_) { + const auto seconds = other.totals_.at(name); + auto it = totals_.find(name); + if (it == totals_.end()) { + order_.push_back(name); + totals_.emplace(name, seconds); + } else { + it->second += seconds; + } + } + } + + void report(const char *heading = "[bioimage profile]") const { + std::fprintf(stderr, "%s\n", heading); double total = 0.0; for (const auto *name : order_) { total += totals_.at(name); @@ -88,7 +104,11 @@ inline ProfileTimerNull make_profile_timer(NullProfiler &profiler, const char *n } // namespace bioimage_cpp::detail -#define BIOIMAGE_PROFILE_INIT(var) ::bioimage_cpp::detail::Profiler var; +namespace bioimage_cpp::detail { +using ActiveProfiler = Profiler; +} // namespace bioimage_cpp::detail + +#define BIOIMAGE_PROFILE_INIT(var) ::bioimage_cpp::detail::ActiveProfiler var; // Two-level indirection so __LINE__ is expanded before the token paste; without // it every BIOIMAGE_PROFILE_SCOPE in a translation unit would declare the same // identifier `_bp___LINE__`, breaking two scopes in one block. @@ -96,11 +116,17 @@ inline ProfileTimerNull make_profile_timer(NullProfiler &profiler, const char *n #define BIOIMAGE_PROFILE_CONCAT(a, b) BIOIMAGE_PROFILE_CONCAT_(a, b) #define BIOIMAGE_PROFILE_SCOPE(var, name) auto BIOIMAGE_PROFILE_CONCAT(_bp_, __LINE__) = ::bioimage_cpp::detail::make_profile_timer(var, name); #define BIOIMAGE_PROFILE_REPORT(var) (var).report(); +#define BIOIMAGE_PROFILE_REPORT_NAMED(var, heading) (var).report(heading); #else -#define BIOIMAGE_PROFILE_INIT(var) ::bioimage_cpp::detail::NullProfiler var; +namespace bioimage_cpp::detail { +using ActiveProfiler = NullProfiler; +} // namespace bioimage_cpp::detail + +#define BIOIMAGE_PROFILE_INIT(var) ::bioimage_cpp::detail::ActiveProfiler var; #define BIOIMAGE_PROFILE_SCOPE(var, name) (void)var; #define BIOIMAGE_PROFILE_REPORT(var) (void)var; +#define BIOIMAGE_PROFILE_REPORT_NAMED(var, heading) (void)var; #endif diff --git a/include/bioimage_cpp/distance/distance_transform.hxx b/include/bioimage_cpp/distance/distance_transform.hxx index e7738ec..9d5e71d 100644 --- a/include/bioimage_cpp/distance/distance_transform.hxx +++ b/include/bioimage_cpp/distance/distance_transform.hxx @@ -20,6 +20,10 @@ namespace detail { constexpr double kInfinity = std::numeric_limits::infinity(); +using SquaredDistanceConsumer = void (*)( + const double *, std::size_t, void * +); + inline std::ptrdiff_t number_of_elements(const std::vector &shape) { std::ptrdiff_t n = 1; for (const auto axis_size : shape) { @@ -186,11 +190,14 @@ struct DistanceTransformOutputs { // >=1 = explicit thread count. Threading splits the orthogonal lines of each // axis sweep across threads via detail::parallel_for_chunks; the per-axis sweep // is a barrier (the next axis depends on the current axis result). -inline void distance_transform( +inline void detail_distance_transform_impl( const ConstArrayView &input, const std::vector &sampling, const DistanceTransformOutputs &outputs, - std::size_t n_threads = 1 + const std::size_t n_threads, + std::unique_ptr initialized_squared_distance, + const detail::SquaredDistanceConsumer squared_distance_consumer, + void *squared_distance_consumer_context ) { const auto ndim = input.ndim(); if (ndim < 1) { @@ -220,6 +227,8 @@ inline void distance_transform( const bool want_distances = outputs.distances.data != nullptr; const bool want_indices = outputs.indices.data != nullptr; const bool want_vectors = outputs.vectors.data != nullptr; + const bool use_initialized_squared = + initialized_squared_distance != nullptr; const bool track_feature = want_indices || want_vectors; bool is_isotropic = true; for (std::size_t axis = 0; axis < sampling.size(); ++axis) { @@ -235,8 +244,8 @@ inline void distance_transform( // and indices against a virtual background row at axis-0 coordinate -1; we // mirror that convention so callers can switch between SciPy and us // without surprises. - bool has_background = false; - { + bool has_background = use_initialized_squared; + if (!use_initialized_squared) { BIOIMAGE_PROFILE_SCOPE(profiler, "scan_for_bg") for (std::ptrdiff_t i = 0; i < n; ++i) { if (input.data[i] == 0) { @@ -277,9 +286,9 @@ inline void distance_transform( // Squared sampled distance buffer. ndim per-axis feature-coord buffers // (int32) replace the previous flat int64 feature index — this lets the // output pass materialize indices/vectors without re-unraveling per pixel. - auto squared_distance = std::make_unique_for_overwrite( - static_cast(n) - ); + auto squared_distance = use_initialized_squared + ? std::move(initialized_squared_distance) + : std::make_unique_for_overwrite(static_cast(n)); std::vector> feature_coord; if (track_feature) { feature_coord.resize(static_cast(ndim)); @@ -315,7 +324,7 @@ inline void distance_transform( // The first-axis gather initializes the uninitialized squared- // distance buffer directly from the input. Later axes gather the // preceding sweep, avoiding a redundant full-volume init pass. - if (ax == 0) { + if (ax == 0 && !use_initialized_squared) { for (std::ptrdiff_t i = 0; i < line_length; ++i) { const auto index = static_cast(base + i * stride); ws.f[static_cast(i)] = @@ -462,7 +471,49 @@ inline void distance_transform( } } } + if (squared_distance_consumer != nullptr) { + BIOIMAGE_PROFILE_SCOPE(profiler, "squared_distance_consumer") + squared_distance_consumer( + squared_distance.get(), static_cast(n), + squared_distance_consumer_context + ); + } BIOIMAGE_PROFILE_REPORT(profiler) } +inline void distance_transform( + const ConstArrayView &input, + const std::vector &sampling, + const DistanceTransformOutputs &outputs, + const std::size_t n_threads = 1 +) { + detail_distance_transform_impl( + input, sampling, outputs, n_threads, nullptr, nullptr, nullptr + ); +} + +namespace detail { + +inline void distance_transform_from_squared( + const std::vector &shape, + const std::vector &sampling, + std::unique_ptr initialized_squared_distance, + const SquaredDistanceConsumer consumer, + void *consumer_context, + const std::size_t n_threads +) { + if (initialized_squared_distance == nullptr) { + throw std::invalid_argument( + "initialized squared-distance buffer must not be null" + ); + } + ConstArrayView shape_only_input{nullptr, shape, {}}; + detail_distance_transform_impl( + shape_only_input, sampling, {}, n_threads, + std::move(initialized_squared_distance), consumer, consumer_context + ); +} + +} // namespace detail + } // namespace bioimage_cpp::distance diff --git a/include/bioimage_cpp/skeleton/detail/compact_grid_dijkstra.hxx b/include/bioimage_cpp/skeleton/detail/compact_grid_dijkstra.hxx index 7806b05..993b2c1 100644 --- a/include/bioimage_cpp/skeleton/detail/compact_grid_dijkstra.hxx +++ b/include/bioimage_cpp/skeleton/detail/compact_grid_dijkstra.hxx @@ -26,6 +26,7 @@ enum class CompactAdjacency { struct CompactNeighbor { std::ptrdiff_t delta = 0; double physical_length = 0.0; + std::array coordinate_delta{}; }; // A deterministic foreground-only view of a zero-padded 3D mask. Compact IDs @@ -94,6 +95,11 @@ inline void build_compact_neighbors( static_cast(dz) * domain.strides[0] + static_cast(dy) * domain.strides[1] + dx, std::sqrt(pz * pz + py * py + px * px), + { + static_cast(dz), + static_cast(dy), + static_cast(dx), + }, }; } } @@ -281,6 +287,31 @@ inline void for_each_compact_neighbor( } } +template +inline void for_each_compact_neighbor_with_metadata( + const CompactGridDomain &domain, + const std::uint32_t node, + const Body &body +) { + if constexpr (Adjacency == CompactAdjacency::Csr) { + for (auto edge = domain.offsets[node]; edge < domain.offsets[node + 1]; ++edge) { + body( + domain.targets[edge], + domain.neighbors[domain.neighbor_codes[edge]] + ); + } + } else { + const auto full = static_cast(domain.compact_to_full[node]); + for (const auto &neighbor : domain.neighbors) { + const auto target_full = static_cast(full + neighbor.delta); + const auto target = domain.full_to_compact[target_full]; + if (target != kNoCompactNode) { + body(target, neighbor); + } + } + } +} + template inline void compact_physical_distance_field( const CompactGridDomain &domain, diff --git a/include/bioimage_cpp/skeleton/detail/components.hxx b/include/bioimage_cpp/skeleton/detail/components.hxx index 8794663..21bf205 100644 --- a/include/bioimage_cpp/skeleton/detail/components.hxx +++ b/include/bioimage_cpp/skeleton/detail/components.hxx @@ -93,6 +93,9 @@ struct PreparedTeasarComponent { // extends foreground across artificial processing-block cuts so those // cuts do not become false object boundaries. std::vector distance_mask; + // Optional foreground-only EDT values in ascending padded C-order. + // Binary dispatch can compute these once for several components. + std::vector compact_dbf; std::array input_origin{}; std::size_t foreground_count = 0; std::vector required_target_voxels; diff --git a/include/bioimage_cpp/skeleton/detail/invalidation.hxx b/include/bioimage_cpp/skeleton/detail/invalidation.hxx index 14fe0b8..5441f1a 100644 --- a/include/bioimage_cpp/skeleton/detail/invalidation.hxx +++ b/include/bioimage_cpp/skeleton/detail/invalidation.hxx @@ -1,11 +1,14 @@ #pragma once #include "bioimage_cpp/detail/grid.hxx" +#include "bioimage_cpp/skeleton/detail/compact_grid_dijkstra.hxx" +#include #include #include #include #include +#include #include #include #include @@ -34,6 +37,231 @@ struct BallInvalidationGreater { } }; +struct CompactBallInvalidationEntry { + double distance = 0.0; + std::uint32_t source = 0; + std::uint32_t node = 0; +}; + +struct CompactBallInvalidationGreater { + bool operator()( + const CompactBallInvalidationEntry &first, + const CompactBallInvalidationEntry &second + ) const noexcept { + if (first.distance != second.distance) { + return first.distance > second.distance; + } + if (first.source != second.source) { + return first.source > second.source; + } + return first.node > second.node; + } +}; + +struct CompactBallInvalidationStats { + std::size_t offers = 0; + std::size_t accepted = 0; + std::size_t pops = 0; + std::size_t stale_pops = 0; + std::size_t invalidated = 0; + std::size_t peak_heap = 0; + + void reset() noexcept { + *this = {}; + } + + void merge(const CompactBallInvalidationStats &other) noexcept { + offers += other.offers; + accepted += other.accepted; + pops += other.pops; + stale_pops += other.stale_pops; + invalidated += other.invalidated; + peak_heap = std::max(peak_heap, other.peak_heap); + } +}; + +struct CompactBallInvalidationWorkspace { + std::vector pending_distance; + std::vector pending_source; + std::vector pending_generation; + std::vector heap; + std::vector> source_coordinates; + std::uint32_t generation = 0; + + void begin(const std::size_t number_of_nodes, const std::size_t number_of_sources) { + if (pending_generation.size() != number_of_nodes) { + pending_distance.resize(number_of_nodes); + pending_source.resize(number_of_nodes); + pending_generation.assign(number_of_nodes, 0); + generation = 0; + } + if (generation == std::numeric_limits::max()) { + std::fill(pending_generation.begin(), pending_generation.end(), 0); + generation = 1; + } else { + ++generation; + } + heap.clear(); + source_coordinates.resize(number_of_sources); + } +}; + +inline void compact_ball_heap_push( + CompactBallInvalidationWorkspace &workspace, + const CompactBallInvalidationEntry entry, + CompactBallInvalidationStats *stats +) { + workspace.heap.push_back(entry); + std::push_heap( + workspace.heap.begin(), workspace.heap.end(), + CompactBallInvalidationGreater{} + ); + if (stats != nullptr) { + ++stats->accepted; + stats->peak_heap = std::max(stats->peak_heap, workspace.heap.size()); + } +} + +inline CompactBallInvalidationEntry compact_ball_heap_pop( + CompactBallInvalidationWorkspace &workspace, + CompactBallInvalidationStats *stats +) { + std::pop_heap( + workspace.heap.begin(), workspace.heap.end(), + CompactBallInvalidationGreater{} + ); + const auto entry = workspace.heap.back(); + workspace.heap.pop_back(); + if (stats != nullptr) { + ++stats->pops; + } + return entry; +} + +// Invalidate strict physical-radius balls in compact foreground space. Keep +// only the best discovered entry per active node. A displaced heap entry is +// stale and cannot affect traversal because the better pending entry sorts +// before it. +template +inline std::size_t invalidate_compact_path_balls( + std::vector &active, + const std::span path, + const std::span radii, + const CompactGridDomain &domain, + const std::array &spacing, + CompactBallInvalidationWorkspace &workspace, + CompactBallInvalidationStats *stats = nullptr +) { + if (domain.shape.size() != 3) { + throw std::invalid_argument("ball invalidation requires a 3D domain"); + } + if (path.size() != radii.size()) { + throw std::invalid_argument("ball invalidation path and radii must match"); + } + if (active.size() != domain.size()) { + throw std::invalid_argument("ball invalidation active domain mismatch"); + } + if constexpr (Adjacency == CompactAdjacency::Csr) { + if (!domain.has_csr()) { + throw std::invalid_argument("compact CSR adjacency is not available"); + } + } else if (!domain.has_full_lookup()) { + throw std::invalid_argument("compact full-index lookup is not available"); + } + if (path.size() > static_cast(kNoCompactNode)) { + throw std::invalid_argument("ball invalidation path exceeds uint32 range"); + } + + workspace.begin(domain.size(), path.size()); + for (std::size_t source = 0; source < path.size(); ++source) { + if (path[source] >= domain.size()) { + throw std::invalid_argument("ball invalidation path node is out of bounds"); + } + bioimage_cpp::detail::coords_from_index( + domain.compact_to_full[path[source]], domain.strides, 3, + workspace.source_coordinates[source].data() + ); + } + + const auto offer = [&](const std::uint32_t node, + const double distance, + const std::uint32_t source) { + if (stats != nullptr) { + ++stats->offers; + } + if (active[node] == 0) { + return; + } + const bool pending = + workspace.pending_generation[node] == workspace.generation; + if ( + pending && + !(distance < workspace.pending_distance[node] || + (distance == workspace.pending_distance[node] && + source < workspace.pending_source[node])) + ) { + return; + } + workspace.pending_distance[node] = distance; + workspace.pending_source[node] = source; + workspace.pending_generation[node] = workspace.generation; + compact_ball_heap_push(workspace, {distance, source, node}, stats); + }; + + for (std::uint32_t source = 0; source < path.size(); ++source) { + offer(path[source], 0.0, source); + } + + std::size_t invalidated = 0; + std::array coordinate{}; + while (!workspace.heap.empty()) { + const auto entry = compact_ball_heap_pop(workspace, stats); + const bool current = + workspace.pending_generation[entry.node] == workspace.generation && + workspace.pending_distance[entry.node] == entry.distance && + workspace.pending_source[entry.node] == entry.source; + if (!current || active[entry.node] == 0) { + if (stats != nullptr) { + ++stats->stale_pops; + } + continue; + } + + active[entry.node] = 0; + ++invalidated; + if (stats != nullptr) { + ++stats->invalidated; + } + bioimage_cpp::detail::coords_from_index( + domain.compact_to_full[entry.node], domain.strides, 3, + coordinate.data() + ); + for_each_compact_neighbor_with_metadata( + domain, entry.node, + [&](const std::uint32_t target, const CompactNeighbor &neighbor) { + if (active[target] == 0) { + return; + } + double distance_squared = 0.0; + for (std::size_t axis = 0; axis < 3; ++axis) { + const auto neighbor_coordinate = + coordinate[axis] + neighbor.coordinate_delta[axis]; + const auto delta = static_cast( + neighbor_coordinate - + workspace.source_coordinates[entry.source][axis] + ) * spacing[axis]; + distance_squared += delta * delta; + } + const double distance = std::sqrt(distance_squared); + if (distance < radii[entry.source]) { + offer(target, distance, entry.source); + } + } + ); + } + return invalidated; +} + // Invalidate the active, 26-connected part of each strict physical-radius ball. // Inactive voxels are barriers. Active source voxels are always invalidated. inline std::size_t invalidate_path_balls( diff --git a/include/bioimage_cpp/skeleton/teasar.hxx b/include/bioimage_cpp/skeleton/teasar.hxx index d8e9ae6..33b3085 100644 --- a/include/bioimage_cpp/skeleton/teasar.hxx +++ b/include/bioimage_cpp/skeleton/teasar.hxx @@ -546,7 +546,9 @@ inline LatticeSkeletonGraph teasar_compact_impl( const ConstArrayView &mask, detail::PreparedTeasarComponent *prepared, const TeasarOptions &options, - const bool report_profile + const bool report_profile, + detail::CompactBallInvalidationStats *ball_invalidation_stats = nullptr, + bioimage_cpp::detail::ActiveProfiler *component_profile = nullptr ) { if (prepared == nullptr) { detail_teasar::validate_options(mask, options); @@ -562,6 +564,7 @@ inline LatticeSkeletonGraph teasar_compact_impl( std::size_t n = 0; std::vector padded_mask; std::vector distance_mask; + std::vector compact_dbf; std::vector required_targets; std::size_t required_root = std::numeric_limits::max(); if (prepared != nullptr) { @@ -570,6 +573,7 @@ inline LatticeSkeletonGraph teasar_compact_impl( shape = std::move(prepared->padded_shape); padded_mask = std::move(prepared->padded_mask); distance_mask = std::move(prepared->distance_mask); + compact_dbf = std::move(prepared->compact_dbf); required_targets = std::move(prepared->required_target_voxels); required_root = prepared->required_root_voxel; n = padded_mask.size(); @@ -634,21 +638,24 @@ inline LatticeSkeletonGraph teasar_compact_impl( options.number_of_threads, foreground_count ); - auto dbf = std::make_unique_for_overwrite(n); - { - BIOIMAGE_PROFILE_SCOPE(profile, "distance_transform") - const auto &distance_input = distance_mask.empty() - ? padded_mask : distance_mask; - ConstArrayView padded_view{ - distance_input.data(), shape, {} - }; - ArrayView distances_view{dbf.get(), shape, {}}; - distance::distance_transform( - padded_view, - {options.spacing[0], options.spacing[1], options.spacing[2]}, - {distances_view, {}, {}}, - effective_threads - ); + std::unique_ptr dbf; + if (compact_dbf.empty()) { + dbf = std::make_unique_for_overwrite(n); + { + BIOIMAGE_PROFILE_SCOPE(profile, "distance_transform") + const auto &distance_input = distance_mask.empty() + ? padded_mask : distance_mask; + ConstArrayView padded_view{ + distance_input.data(), shape, {} + }; + ArrayView distances_view{dbf.get(), shape, {}}; + distance::distance_transform( + padded_view, + {options.spacing[0], options.spacing[1], options.spacing[2]}, + {distances_view, {}, {}}, + effective_threads + ); + } } detail::CompactGridDomain domain; @@ -678,16 +685,26 @@ inline LatticeSkeletonGraph teasar_compact_impl( } double dbf_max = 0.0; - std::vector compact_dbf; - compact_dbf.reserve(domain.size()); - { - BIOIMAGE_PROFILE_SCOPE(profile, "dbf_compaction") - for (std::uint32_t node = 0; node < domain.size(); ++node) { - const auto value = dbf[domain.compact_to_full[node]]; - compact_dbf.push_back(value); + if (compact_dbf.empty()) { + compact_dbf.reserve(domain.size()); + { + BIOIMAGE_PROFILE_SCOPE(profile, "dbf_compaction") + for (std::uint32_t node = 0; node < domain.size(); ++node) { + const auto value = dbf[domain.compact_to_full[node]]; + compact_dbf.push_back(value); + dbf_max = std::max(dbf_max, static_cast(value)); + } + dbf.reset(); + } + } else { + if (compact_dbf.size() != domain.size()) { + throw std::runtime_error( + "TEASAR precomputed distance count is inconsistent" + ); + } + for (const auto value : compact_dbf) { dbf_max = std::max(dbf_max, static_cast(value)); } - dbf.reset(); } detail::CompactDijkstraWorkspace dijkstra_workspace; @@ -759,7 +776,13 @@ inline LatticeSkeletonGraph teasar_compact_impl( } } - std::vector active = std::move(padded_mask); + std::vector active; + if (options.invalidation == TeasarInvalidation::Ball) { + active.assign(domain.size(), std::uint8_t{1}); + std::vector().swap(padded_mask); + } else { + active = std::move(padded_mask); + } std::size_t active_count = foreground_count; std::vector vertex_of_node(domain.size(), -1); std::vector skeleton_nodes; @@ -795,6 +818,8 @@ inline LatticeSkeletonGraph teasar_compact_impl( add_vertex(root); std::vector path; + std::vector ball_radii; + detail::CompactBallInvalidationWorkspace ball_invalidation_workspace; detail::RowIntervalUnion invalidated_rows( n / static_cast(shape[2]), shape[2] ); @@ -836,24 +861,20 @@ inline LatticeSkeletonGraph teasar_compact_impl( { BIOIMAGE_PROFILE_SCOPE(profile, "invalidation") if (options.invalidation == TeasarInvalidation::Ball) { - std::vector full_path; - std::vector radii; - full_path.reserve(path.size()); - radii.reserve(path.size()); + ball_radii.clear(); + ball_radii.reserve(path.size()); for (const auto node : path) { - full_path.push_back( - static_cast(domain.compact_to_full[node]) - ); const double radius = options.scale * static_cast(compact_dbf[node]) + options.constant; if (!std::isfinite(radius)) { throw std::runtime_error("TEASAR invalidation radius overflowed"); } - radii.push_back(radius); + ball_radii.push_back(radius); } - const auto invalidated = detail::invalidate_path_balls( - active, full_path, radii, shape, options.spacing + const auto invalidated = detail::invalidate_compact_path_balls( + active, path, ball_radii, domain, options.spacing, + ball_invalidation_workspace, ball_invalidation_stats ); if (invalidated > active_count) { throw std::runtime_error( @@ -963,8 +984,11 @@ inline LatticeSkeletonGraph teasar_compact_impl( ) { Distance target_distance = Distance{-1}; for (std::uint32_t node = 0; node < domain.size(); ++node) { - const auto full = domain.compact_to_full[node]; - if (active[full] != 0 && root_field[node] > target_distance) { + const auto active_index = options.invalidation == + TeasarInvalidation::Ball + ? static_cast(node) + : static_cast(domain.compact_to_full[node]); + if (active[active_index] != 0 && root_field[node] > target_distance) { target = node; target_distance = root_field[node]; } @@ -974,7 +998,11 @@ inline LatticeSkeletonGraph teasar_compact_impl( if (!targets_ordered) { ordered_targets.reserve(active_count); for (std::uint32_t node = 0; node < domain.size(); ++node) { - if (active[domain.compact_to_full[node]] != 0) { + const auto active_index = options.invalidation == + TeasarInvalidation::Ball + ? static_cast(node) + : static_cast(domain.compact_to_full[node]); + if (active[active_index] != 0) { ordered_targets.push_back(node); } } @@ -992,7 +1020,13 @@ inline LatticeSkeletonGraph teasar_compact_impl( while ( ordered_target_cursor < ordered_targets.size() && active[ - domain.compact_to_full[ordered_targets[ordered_target_cursor]] + options.invalidation == TeasarInvalidation::Ball + ? static_cast( + ordered_targets[ordered_target_cursor] + ) + : static_cast(domain.compact_to_full[ + ordered_targets[ordered_target_cursor] + ]) ] == 0 ) { ++ordered_target_cursor; @@ -1011,6 +1045,9 @@ inline LatticeSkeletonGraph teasar_compact_impl( if (report_profile) { BIOIMAGE_PROFILE_REPORT(profile) } + if (component_profile != nullptr) { + component_profile->merge(profile); + } return graph; } @@ -1027,11 +1064,14 @@ inline LatticeSkeletonGraph teasar_compact( template inline LatticeSkeletonGraph teasar_compact_prepared( detail::PreparedTeasarComponent prepared, - const TeasarOptions &options + const TeasarOptions &options, + detail::CompactBallInvalidationStats *ball_invalidation_stats = nullptr, + bioimage_cpp::detail::ActiveProfiler *component_profile = nullptr ) { const ConstArrayView unused{}; return teasar_compact_impl( - unused, &prepared, options, false + unused, &prepared, options, false, ball_invalidation_stats, + component_profile ); } @@ -1235,19 +1275,341 @@ std::vector component_thread_budgets( return budgets; } +enum class ComponentEdtStrategy { + Auto, + Local, + Shared, +}; + +enum class SharedEdtDecision { + Selected, + ForcedLocal, + FewerThanTwoComponents, + CompactRange, + VolumeRatio, + ScratchLimit, +}; + +inline constexpr std::size_t kSharedEdtScratchLimitBytes = + std::size_t{256} * 1024 * 1024; + +inline const char *shared_edt_decision_name(const SharedEdtDecision decision) { + switch (decision) { + case SharedEdtDecision::Selected: + return "shared"; + case SharedEdtDecision::ForcedLocal: + return "forced-local"; + case SharedEdtDecision::FewerThanTwoComponents: + return "fewer-than-two-components"; + case SharedEdtDecision::CompactRange: + return "compact-range"; + case SharedEdtDecision::VolumeRatio: + return "volume-ratio"; + case SharedEdtDecision::ScratchLimit: + return "scratch-limit"; + } + return "unknown"; +} + +struct SharedEdtPreparation { + std::vector> component_dbf; + SharedEdtDecision decision = SharedEdtDecision::ForcedLocal; + std::size_t estimated_scratch_bytes = 0; + + [[nodiscard]] bool selected() const noexcept { + return decision == SharedEdtDecision::Selected; + } +}; + +struct SharedEdtGatherContext { + const detail::ComponentSet *components = nullptr; + std::array global_begin{}; + const std::vector *shared_shape = nullptr; + std::vector> *component_dbf = nullptr; +}; + +inline void gather_shared_squared_distances( + const double *squared_distances, + const std::size_t number_of_values, + void *raw_context +) { + auto &context = *static_cast(raw_context); + if ( + context.components == nullptr || context.shared_shape == nullptr || + context.component_dbf == nullptr + ) { + throw std::invalid_argument("shared EDT gather context is incomplete"); + } + const auto expected_values = detail::checked_shape_size( + *context.shared_shape, "shared EDT gather shape overflows size_t" + ); + if (number_of_values != expected_values) { + throw std::runtime_error("shared EDT gather volume is inconsistent"); + } + const auto strides = bioimage_cpp::detail::c_order_strides( + *context.shared_shape + ); + const auto &components = *context.components; + auto &component_dbf = *context.component_dbf; + for (std::size_t component_id = 0; + component_id < components.components.size(); ++component_id) { + const auto &component = components.components[component_id]; + auto &values = component_dbf[component_id]; + values.reserve(static_cast(component.voxel_count)); + for (std::size_t offset = 0; offset < component.number_of_runs; ++offset) { + const auto run_id = components.component_run_ids[ + component.run_offset + offset + ]; + const auto &run = components.runs[run_id]; + const auto z = run.z - context.global_begin[0] + 1; + const auto y = run.y - context.global_begin[1] + 1; + const auto row_begin = static_cast( + z * strides[0] + y * strides[1] + ); + for (auto x = run.x_begin; x <= run.x_end; ++x) { + const auto local_x = static_cast( + x - context.global_begin[2] + 1 + ); + values.push_back(static_cast( + std::sqrt(squared_distances[row_begin + local_x]) + )); + } + } + if (values.size() != static_cast(component.voxel_count)) { + throw std::runtime_error( + "shared EDT component value count is inconsistent" + ); + } + } +} + +inline std::size_t shared_edt_scratch_estimate( + const std::size_t shared_volume, + const std::size_t foreground_count, + const std::size_t maximum_extent, + const std::size_t number_of_threads +) { + const auto volume_edt_bytes = detail::checked_multiply_size( + shared_volume, sizeof(double), + "shared EDT volume scratch estimate overflows size_t" + ); + const auto compact_output_bytes = detail::checked_multiply_size( + foreground_count, sizeof(float), + "shared EDT compact output estimate overflows size_t" + ); + const auto workspace_line_bytes = detail::checked_add_size( + detail::checked_multiply_size( + maximum_extent, std::size_t{32}, + "shared EDT line scratch estimate overflows size_t" + ), + std::size_t{8 + 6 * sizeof(std::vector)}, + "shared EDT line scratch estimate overflows size_t" + ); + const auto workspace_bytes = detail::checked_multiply_size( + workspace_line_bytes, number_of_threads, + "shared EDT thread scratch estimate overflows size_t" + ); + const auto edt_peak = detail::checked_add_size( + detail::checked_add_size( + volume_edt_bytes, compact_output_bytes, + "shared EDT peak estimate overflows size_t" + ), + workspace_bytes, + "shared EDT peak estimate overflows size_t" + ); + const auto compact_dbf_bytes = detail::checked_multiply_size( + foreground_count, sizeof(float), + "shared EDT compact DBF estimate overflows size_t" + ); + const auto gather_peak = detail::checked_multiply_size( + compact_dbf_bytes, std::size_t{1}, + "shared EDT gather peak estimate overflows size_t" + ); + const auto raw_peak = std::max(edt_peak, gather_peak); + const auto headroom = detail::checked_add_size( + raw_peak, std::size_t{9}, + "shared EDT headroom estimate overflows size_t" + ) / 10; + return detail::checked_add_size( + raw_peak, headroom, + "shared EDT scratch estimate overflows size_t" + ); +} + +template +SharedEdtPreparation prepare_shared_binary_edt( + const detail::ComponentSet &components, + const TeasarOptions &options, + const ComponentEdtStrategy strategy, + const std::size_t scratch_limit_bytes, + Profiler &profile +) { + SharedEdtPreparation result; + const auto count = components.components.size(); + if (strategy == ComponentEdtStrategy::Local) { + result.decision = SharedEdtDecision::ForcedLocal; + return result; + } + if (count == 0) { + result.decision = SharedEdtDecision::FewerThanTwoComponents; + return result; + } + if (strategy == ComponentEdtStrategy::Auto && count < 2) { + result.decision = SharedEdtDecision::FewerThanTwoComponents; + return result; + } + + auto global_begin = components.components.front().begin; + auto global_end = components.components.front().end; + std::size_t local_volume_sum = 0; + for (const auto &component : components.components) { + const auto padded_volume = detail::padded_component_volume(component); + if ( + padded_volume > static_cast(detail::kNoCompactNode) || + component.voxel_count > detail::kNoCompactNode + ) { + if (strategy == ComponentEdtStrategy::Shared) { + throw std::invalid_argument( + "forced shared EDT requires uint32-compatible components" + ); + } + result.decision = SharedEdtDecision::CompactRange; + return result; + } + local_volume_sum = detail::checked_add_size( + local_volume_sum, padded_volume, + "summed component EDT volume overflows size_t" + ); + for (std::size_t axis = 0; axis < 3; ++axis) { + global_begin[axis] = std::min(global_begin[axis], component.begin[axis]); + global_end[axis] = std::max(global_end[axis], component.end[axis]); + } + } + + std::vector shared_shape(3); + std::size_t maximum_extent = 0; + for (std::size_t axis = 0; axis < 3; ++axis) { + const auto extent = global_end[axis] - global_begin[axis]; + if (extent > std::numeric_limits::max() - 2) { + throw std::overflow_error("shared EDT shape overflows ptrdiff_t"); + } + shared_shape[axis] = extent + 2; + maximum_extent = std::max( + maximum_extent, static_cast(shared_shape[axis]) + ); + } + const auto shared_volume = detail::checked_shape_size( + shared_shape, "shared EDT volume overflows size_t" + ); + if (shared_volume > static_cast(detail::kNoCompactNode)) { + if (strategy == ComponentEdtStrategy::Shared) { + throw std::invalid_argument( + "forced shared EDT volume exceeds uint32 index range" + ); + } + result.decision = SharedEdtDecision::CompactRange; + return result; + } + if ( + strategy == ComponentEdtStrategy::Auto && + static_cast(shared_volume) > + 0.75L * static_cast(local_volume_sum) + ) { + result.decision = SharedEdtDecision::VolumeRatio; + return result; + } + + const auto effective_threads = bioimage_cpp::detail::normalize_thread_count( + options.number_of_threads, components.foreground_count + ); + result.estimated_scratch_bytes = shared_edt_scratch_estimate( + shared_volume, components.foreground_count, maximum_extent, + effective_threads + ); + if ( + strategy == ComponentEdtStrategy::Auto && + result.estimated_scratch_bytes > scratch_limit_bytes + ) { + result.decision = SharedEdtDecision::ScratchLimit; + return result; + } + + std::unique_ptr initialized_squared; + { + BIOIMAGE_PROFILE_SCOPE(profile, "shared_edt_setup") + initialized_squared = std::make_unique_for_overwrite( + shared_volume + ); + std::fill( + initialized_squared.get(), + initialized_squared.get() + shared_volume, + 0.0 + ); + const auto strides = bioimage_cpp::detail::c_order_strides(shared_shape); + for (const auto &run : components.runs) { + const auto z = run.z - global_begin[0] + 1; + const auto y = run.y - global_begin[1] + 1; + const auto x_begin = run.x_begin - global_begin[2] + 1; + const auto x_end = run.x_end - global_begin[2] + 1; + const auto row_begin = static_cast( + z * strides[0] + y * strides[1] + ); + for (auto x = x_begin; x <= x_end; ++x) { + const auto index = row_begin + static_cast(x); + initialized_squared[index] = distance::detail::kInfinity; + } + } + } + + result.component_dbf.resize(count); + SharedEdtGatherContext gather_context{ + &components, global_begin, &shared_shape, &result.component_dbf + }; + { + BIOIMAGE_PROFILE_SCOPE(profile, "shared_distance_transform") + distance::detail::distance_transform_from_squared( + shared_shape, + {options.spacing[0], options.spacing[1], options.spacing[2]}, + std::move(initialized_squared), + &gather_shared_squared_distances, &gather_context, + effective_threads + ); + } + result.decision = SharedEdtDecision::Selected; + return result; +} + template std::vector skeletonize_components( const detail::ComponentSet &components, const TeasarOptions &options, const bool include_label_in_errors, const std::vector>> *required_targets = nullptr, - const detail::OpenBlockFaces *open_faces = nullptr + const detail::OpenBlockFaces *open_faces = nullptr, + std::vector> *precomputed_component_dbf = nullptr, + std::vector *ball_invalidation_stats = nullptr, + bioimage_cpp::detail::ActiveProfiler *component_profile = nullptr ) { const auto count = components.components.size(); std::vector results(count); if (count == 0) { return results; } + if (ball_invalidation_stats != nullptr) { + ball_invalidation_stats->assign(count, {}); + } + std::vector component_profiles; + if (component_profile != nullptr) { + component_profiles.resize(count); + } + const auto merge_component_profiles = [&] { + if (component_profile == nullptr) { + return; + } + for (const auto &local_profile : component_profiles) { + component_profile->merge(local_profile); + } + }; const auto total_budget = bioimage_cpp::detail::normalize_thread_count( options.number_of_threads, components.foreground_count ); @@ -1281,11 +1643,22 @@ std::vector skeletonize_components( components, component_id, required_targets->at(component_id), open_faces ); + if (precomputed_component_dbf != nullptr) { + prepared.compact_dbf = std::move( + precomputed_component_dbf->at(component_id) + ); + } auto local_options = options; local_options.number_of_threads = local_budget; results[component_id] = teasar_compact_prepared< detail::CompactAdjacency::OnTheFly, double - >(std::move(prepared), local_options); + >( + std::move(prepared), local_options, + ball_invalidation_stats == nullptr + ? nullptr : &ball_invalidation_stats->at(component_id), + component_profile == nullptr + ? nullptr : &component_profiles[component_id] + ); } catch (const std::exception &error) { throw std::runtime_error( component_context( @@ -1297,6 +1670,7 @@ std::vector skeletonize_components( if (count == 1) { run_component(0, total_budget); + merge_component_profiles(); return results; } if (count >= total_budget) { @@ -1313,6 +1687,7 @@ std::vector skeletonize_components( } } ); + merge_component_profiles(); return results; } @@ -1326,6 +1701,7 @@ std::vector skeletonize_components( } } ); + merge_component_profiles(); return results; } @@ -1360,18 +1736,44 @@ inline SkeletonGraph teasar_with_backend( ); } -inline SkeletonGraph teasar( +inline SkeletonGraph teasar_with_component_edt( const ConstArrayView &mask, - const TeasarOptions &options = {} + const TeasarOptions &options, + const detail_teasar::ComponentEdtStrategy edt_strategy, + const std::size_t shared_edt_scratch_limit, + detail_teasar::SharedEdtDecision *used_decision = nullptr, + std::size_t *estimated_scratch_bytes = nullptr ) { detail_teasar::validate_options(mask, options); BIOIMAGE_PROFILE_INIT(profile) auto components = detail::extract_binary_components(mask, profile); + auto shared_edt = detail_teasar::prepare_shared_binary_edt( + components, options, edt_strategy, shared_edt_scratch_limit, profile + ); + if (used_decision != nullptr) { + *used_decision = shared_edt.decision; + } + if (estimated_scratch_bytes != nullptr) { + *estimated_scratch_bytes = shared_edt.estimated_scratch_bytes; + } std::vector results; + std::vector ball_stats; + bioimage_cpp::detail::ActiveProfiler component_profile; +#ifdef BIOIMAGE_PROFILE + auto *ball_stats_output = &ball_stats; + auto *component_profile_output = &component_profile; +#else + auto *ball_stats_output = + static_cast *>(nullptr); + auto *component_profile_output = + static_cast(nullptr); +#endif { BIOIMAGE_PROFILE_SCOPE(profile, "component_teasar") results = detail_teasar::skeletonize_components( - components, options, false + components, options, false, nullptr, nullptr, + shared_edt.selected() ? &shared_edt.component_dbf : nullptr, + ball_stats_output, component_profile_output ); } std::vector component_ids(results.size()); @@ -1382,11 +1784,50 @@ inline SkeletonGraph teasar( output = detail_teasar::assemble_skeleton_graphs(results, component_ids); } BIOIMAGE_PROFILE_REPORT(profile) + BIOIMAGE_PROFILE_REPORT_NAMED( + component_profile, "[bioimage TEASAR component profile]" + ) +#ifdef BIOIMAGE_PROFILE + detail::CompactBallInvalidationStats combined_ball_stats; + for (const auto &stats : ball_stats) { + combined_ball_stats.merge(stats); + } + std::fprintf( + stderr, + "[bioimage TEASAR diagnostics]\n" + " shared_edt %s\n" + " shared_scratch %zu bytes\n" + " ball_offers %zu\n" + " ball_accepted %zu\n" + " ball_pops %zu\n" + " ball_stale_pops %zu\n" + " ball_invalidated %zu\n" + " ball_peak_heap %zu\n", + detail_teasar::shared_edt_decision_name(shared_edt.decision), + shared_edt.estimated_scratch_bytes, + combined_ball_stats.offers, + combined_ball_stats.accepted, + combined_ball_stats.pops, + combined_ball_stats.stale_pops, + combined_ball_stats.invalidated, + combined_ball_stats.peak_heap + ); +#endif return detail_teasar::lattice_to_physical( std::move(output), options.spacing ); } +inline SkeletonGraph teasar( + const ConstArrayView &mask, + const TeasarOptions &options = {} +) { + return teasar_with_component_edt( + mask, options, detail_teasar::ComponentEdtStrategy::Auto, + detail_teasar::kSharedEdtScratchLimitBytes + ); +} + inline LatticeSkeletonGraph teasar_block( const ConstArrayView &mask, std::vector required_targets, diff --git a/src/bindings/skeleton.cxx b/src/bindings/skeleton.cxx index 7ed71d0..ad736c8 100644 --- a/src/bindings/skeleton.cxx +++ b/src/bindings/skeleton.cxx @@ -143,6 +143,77 @@ nb::tuple teasar_uint8_backend( ); } +nb::tuple teasar_uint8_edt_backend( + UInt8Input mask, + const std::vector &spacing, + const double scale, + const double constant, + const double pdrf_scale, + const double pdrf_exponent, + const bool ball_invalidation, + const bool fix_branching, + const std::string &edt_strategy, + const std::size_t scratch_limit_bytes, + const std::size_t n_threads +) { + if (mask.ndim() != 3) { + throw std::invalid_argument( + "mask must have ndim 3, got ndim=" + std::to_string(mask.ndim()) + ); + } + if (spacing.size() != 3) { + throw std::invalid_argument( + "spacing must contain exactly three values, got " + + std::to_string(spacing.size()) + ); + } + skeleton::detail_teasar::ComponentEdtStrategy selected; + if (edt_strategy == "auto") { + selected = skeleton::detail_teasar::ComponentEdtStrategy::Auto; + } else if (edt_strategy == "local") { + selected = skeleton::detail_teasar::ComponentEdtStrategy::Local; + } else if (edt_strategy == "shared") { + selected = skeleton::detail_teasar::ComponentEdtStrategy::Shared; + } else { + throw std::invalid_argument( + "unknown TEASAR EDT development strategy: " + edt_strategy + ); + } + + std::vector shape(mask.ndim()); + for (std::size_t axis = 0; axis < mask.ndim(); ++axis) { + shape[axis] = static_cast(mask.shape(axis)); + } + ConstArrayView mask_view{mask.data(), shape, {}}; + skeleton::SkeletonGraph result; + skeleton::detail_teasar::SharedEdtDecision decision; + std::size_t estimated_scratch_bytes = 0; + { + nb::gil_scoped_release release; + const skeleton::TeasarOptions options{ + {spacing[0], spacing[1], spacing[2]}, + scale, + constant, + pdrf_scale, + pdrf_exponent, + n_threads, + ball_invalidation ? skeleton::TeasarInvalidation::Ball + : skeleton::TeasarInvalidation::Cube, + fix_branching, + }; + result = skeleton::teasar_with_component_edt( + mask_view, options, selected, scratch_limit_bytes, + &decision, &estimated_scratch_bytes + ); + } + auto arrays = skeleton_graph_to_tuple(result); + return nb::make_tuple( + arrays[0], arrays[1], arrays[2], + skeleton::detail_teasar::shared_edt_decision_name(decision), + estimated_scratch_bytes + ); +} + template nb::dict teasar_labels_impl( nb::ndarray labels, @@ -236,6 +307,23 @@ void bind_skeleton(nb::module_ &m) { nb::arg("n_threads") = 1, "Development-only TEASAR backend selector." ); + m.def( + "_teasar_uint8_edt_backend", + &teasar_uint8_edt_backend, + nb::arg("mask"), + nb::arg("spacing"), + nb::arg("scale"), + nb::arg("constant"), + nb::arg("pdrf_scale"), + nb::arg("pdrf_exponent"), + nb::arg("ball_invalidation"), + nb::arg("fix_branching"), + nb::arg("edt_strategy"), + nb::arg("scratch_limit_bytes") = + skeleton::detail_teasar::kSharedEdtScratchLimitBytes, + nb::arg("n_threads") = 1, + "Development-only TEASAR component EDT selector." + ); #define BIC_BIND_TEASAR_LABELS(name, type) \ m.def( \ diff --git a/tests/skeleton/test_teasar.py b/tests/skeleton/test_teasar.py index 245c44c..888785f 100644 --- a/tests/skeleton/test_teasar.py +++ b/tests/skeleton/test_teasar.py @@ -28,6 +28,31 @@ def _teasar_backend( ) +def _teasar_edt_backend( + mask, + strategy, + *, + spacing=(1.5, 1.0, 1.0), + ball_invalidation=True, + fix_branching=True, + scratch_limit_bytes=256 * 1024 * 1024, + number_of_threads=1, +): + return _core._teasar_uint8_edt_backend( + np.ascontiguousarray(mask, dtype=np.uint8), + spacing, + 1.5, + 1.0, + 100000.0, + 4.0, + ball_invalidation, + fix_branching, + strategy, + scratch_limit_bytes, + number_of_threads, + ) + + def _assert_valid_tree(mask, vertices, edges, radii, spacing=(1.0, 1.0, 1.0)): mask = np.asarray(mask) != 0 spacing = np.asarray(spacing, dtype=np.float64) @@ -296,6 +321,92 @@ def test_compact_fp64_backends_have_exact_dense_parity( np.testing.assert_array_equal(got, expected) +@pytest.mark.parametrize("seed", [7, 19, 41]) +@pytest.mark.parametrize("spacing", [(1.0, 1.0, 1.0), (2.5, 1.25, 0.75)]) +@pytest.mark.parametrize("fix_branching", [True, False]) +def test_coalesced_ball_queue_preserves_randomized_dense_parity( + seed, spacing, fix_branching +): + rng = np.random.default_rng(seed) + mask = np.zeros((17, 19, 21), dtype=np.uint8) + coordinate = np.array([8, 9, 0], dtype=np.int64) + for _ in range(180): + lo = np.maximum(coordinate - 1, 0) + hi = np.minimum(coordinate + 2, mask.shape) + mask[tuple(slice(int(a), int(b)) for a, b in zip(lo, hi))] = 1 + coordinate += rng.integers(-1, 2, size=3) + coordinate = np.clip(coordinate, 0, np.asarray(mask.shape) - 1) + + dense = _teasar_backend( + mask, + "dense-fp64", + spacing=spacing, + ball_invalidation=True, + fix_branching=fix_branching, + ) + for backend in ("compact-on-the-fly-fp64", "compact-csr-fp64"): + compact = _teasar_backend( + mask, + backend, + spacing=spacing, + ball_invalidation=True, + fix_branching=fix_branching, + ) + for got, expected in zip(compact, dense): + np.testing.assert_array_equal(got, expected) + + +def _clustered_components_mask(): + mask = np.zeros((37, 37, 37), dtype=np.uint8) + for lo, hi in ((1, 35), (6, 30), (11, 25), (16, 20)): + mask[lo, lo:hi + 1, lo:hi + 1] = 1 + mask[hi, lo:hi + 1, lo:hi + 1] = 1 + mask[lo:hi + 1, lo, lo:hi + 1] = 1 + mask[lo:hi + 1, hi, lo:hi + 1] = 1 + mask[lo:hi + 1, lo:hi + 1, lo] = 1 + mask[lo:hi + 1, lo:hi + 1, hi] = 1 + return mask + + +@pytest.mark.parametrize("spacing", [(1.0, 1.0, 1.0), (2.5, 1.25, 0.75)]) +@pytest.mark.parametrize("fix_branching", [True, False]) +def test_shared_component_edt_preserves_local_result_exactly( + spacing, fix_branching +): + mask = _clustered_components_mask() + local = _teasar_edt_backend( + mask, "local", spacing=spacing, fix_branching=fix_branching + ) + shared = _teasar_edt_backend( + mask, "shared", spacing=spacing, fix_branching=fix_branching, + number_of_threads=4, + ) + assert local[3] == "forced-local" + assert shared[3] == "shared" + assert shared[4] > 0 + for got, expected in zip(shared[:3], local[:3]): + np.testing.assert_array_equal(got, expected) + + +def test_shared_component_edt_auto_selection_and_memory_fallback(): + mask = _clustered_components_mask() + shared = _teasar_edt_backend(mask, "auto") + limited = _teasar_edt_backend(mask, "auto", scratch_limit_bytes=1) + assert shared[3] == "shared" + assert limited[3] == "scratch-limit" + assert limited[4] == shared[4] + for got, expected in zip(limited[:3], shared[:3]): + np.testing.assert_array_equal(got, expected) + + +def test_shared_component_edt_rejects_unprofitable_union_volume(): + mask = np.zeros((35, 35, 35), dtype=np.uint8) + mask[1, 1, 1:5] = 1 + mask[-2, -2, -5:-1] = 1 + result = _teasar_edt_backend(mask, "auto") + assert result[3] == "volume-ratio" + + @pytest.mark.parametrize("spacing", [(1.0, 1.0, 1.0), (2.5, 1.25, 0.75)]) def test_cropped_compact_backends_preserve_dense_parity(spacing): mask = np.zeros((31, 37, 45), dtype=np.uint8) From ea2008dc6b57eca126e3c0a0d2bc224d4d2e7963 Mon Sep 17 00:00:00 2001 From: Constantin Pape Date: Wed, 5 Aug 2026 08:32:03 +0200 Subject: [PATCH 3/3] Add results from further optimization --- development/skeleton/PERFORMANCE_NOTES.md | 69 ++++++ .../skeleton/TEASAR_OPTIONS_COMPARISON.md | 198 ++++++++++++++++++ .../skeleton/detail/compact_grid_dijkstra.hxx | 70 ++++++- include/bioimage_cpp/skeleton/teasar.hxx | 54 ++++- 4 files changed, 377 insertions(+), 14 deletions(-) create mode 100644 development/skeleton/TEASAR_OPTIONS_COMPARISON.md diff --git a/development/skeleton/PERFORMANCE_NOTES.md b/development/skeleton/PERFORMANCE_NOTES.md index c02ebfe..dd77f26 100644 --- a/development/skeleton/PERFORMANCE_NOTES.md +++ b/development/skeleton/PERFORMANCE_NOTES.md @@ -958,3 +958,72 @@ The example MRC skeleton statistics remained unchanged through both optimization stages. Worker-count checks also remained array-exact. Final verification used the normal build with profile-only instrumentation disabled: `1421 passed`. + +## Priority-queue follow-up (2026-08-04) + +This pass evaluated the next two optimization targets. Neither candidate met +its retention gate, so the production queues remain unchanged. + +### Compact Dijkstra queues + +A standalone C++ probe used the largest component from the centered `200^3` +MRC crop. The padded component shape was `(62, 190, 202)`, with 102,427 +foreground voxels. Each value is the median of nine runs after one warmup. +All candidates produced exact distance and predecessor arrays. + +| queue | physical field | change | parental field | change | +| --- | ---: | ---: | ---: | ---: | +| binary heap | 27.31 ms | reference | 21.76 ms | reference | +| 4-ary heap | 27.02 ms | -1.1% | 22.35 ms | +2.7% | +| radix heap | 26.79 ms | -1.9% | 23.61 ms | +8.5% | + +The 10% primitive-speed gate rejected both candidates before integration. +The physical field settled 102,427 nodes from 109,115 queue entries. Only +6.1% of its entries were stale. The parental field had no stale entries. + +Profile builds now report aggregate compact Dijkstra counters. On ball/fix, +the MRC crop used 344 Dijkstra calls. The calls removed 1,950,359 entries and +rejected 115,819 stale entries. This 5.9% stale-pop rate confirms that queue +replacement has limited headroom in the full workload. + +### Indexed ball-invalidation queue + +The second candidate replaced lazy invalidation entries with +`DenseIndexedHeap`. The candidate kept one entry per active node and updated +its priority in place. It preserved the `(distance, source, node)` order and +passed the focused 101-test TEASAR suite. + +On ball/fix, the indexed queue reduced queue removals from 693,630 to 585,875 +and reduced the peak queue size from 18,059 to 15,022. It replaced 107,755 +stale removals with in-place priority updates. The extra locator writes during +heap swaps offset this reduction. + +The real-mask table compares separate normal builds. Each value is the median +of seven calls after two warmups. + +| setting | workers | lazy queue | indexed queue | change | +| --- | ---: | ---: | ---: | ---: | +| ball/fix | 1 | 974.9 ms | 982.0 ms | +0.7% | +| ball/fix | 8 | 302.3 ms | 289.7 ms | -4.2% | +| ball/parental | 1 | 956.6 ms | 975.8 ms | +2.0% | +| ball/parental | 8 | 285.0 ms | 297.5 ms | +4.4% | + +The one-worker result fails the required 10% improvement for both branching +modes. The eight-worker parental result also regresses by more than 3%. + +The synthetic matrix showed workload-dependent gains. At `256^3`, ball/fix +changed by -2.7% with one worker and -12.5% with eight workers. +Ball/parental changed by -9.1% and -2.0%. The `192^3` ball/parental result +regressed by 4.6% with eight workers. These mixed results do not justify the +additional indexed-heap state. + +Median incremental peak RSS stayed at 52.7 MiB with one worker. The indexed +candidate used 76.3 MiB for ball/fix and 73.9 MiB for ball/parental with eight +workers. These values satisfy the memory gate but do not change the timing +decision. + +The final code retains only the profile-only Dijkstra counters. It removes the +4-ary, radix, and indexed invalidation candidates. The distance-transform +implementation remains outside this optimization pass. Final verification +used the normal build and passed all 1,421 tests. The MRC vertices, edges, +radii, and worker-count comparisons remained array-exact. diff --git a/development/skeleton/TEASAR_OPTIONS_COMPARISON.md b/development/skeleton/TEASAR_OPTIONS_COMPARISON.md new file mode 100644 index 0000000..bf02bf0 --- /dev/null +++ b/development/skeleton/TEASAR_OPTIONS_COMPARISON.md @@ -0,0 +1,198 @@ +# TEASAR option comparison on the example MRC mask + +This report compares the new bioimage-cpp TEASAR settings with Kimimaro 5.8.1. +It uses the example mask and centered crop from the initial skeleton diagnostic. + +## Summary + +- Ball invalidation moves the bioimage-cpp topology closer to Kimimaro. +- `fix_branching=False` reduces degree-3 nodes and short terminal spurs in both + implementations. +- Bioimage-cpp ball modes are 1.9--2.0x faster than Kimimaro with one worker. + They are approximately 3.5x faster with eight workers on this crop. +- Coalesced invalidation and a guarded shared distance transform reduced the + one-worker ball time by 52--53% without changing the measured skeletons. +- The default bioimage-cpp behavior remains cube invalidation with + `fix_branching=True`. + +The outputs do not match exactly. Bioimage-cpp and Kimimaro still differ in +component handling, root selection, distance fields, and path tie-breaking. + +## Measurement setup + +Measured on 2026-08-04. + +The input is `examples/skeleton/00004_gt_mask.mrc`. The source shape is +`(324, 1251, 1251)`. The benchmark uses the centered crop with origin +`(62, 525, 525)`, shape `(200, 200, 200)`, and 585,875 foreground voxels. + +The parameters are: + +| parameter | cube modes | ball modes and Kimimaro | +| --- | ---: | ---: | +| spacing | 10 Å | 10 Å | +| `scale` | 0 | 0 | +| `constant` | 70 Å | 140 Å | +| `pdrf_scale` | 100,000 | 100,000 | +| `pdrf_exponent` | 4 | 4 | + +These constants reproduce the working radii from the initial investigation. +The cube and ball timing comparison therefore includes both a geometry change +and a radius change. The bioimage-cpp ball and Kimimaro rows use identical +TEASAR parameters. + +Kimimaro runs with soma handling disabled. Border fixing, hole filling, dust +filtering, and progress reporting are also disabled. `fix` means +`fix_branching=True`. `parental` means `fix_branching=False`. + +The host has an Intel Core i7-1185G7 with four cores and eight hardware threads. +The environment uses Python 3.13.13, NumPy 2.4.6, bioimage-cpp 0.8.0, +Kimimaro 5.8.1, and EDT 3.1.1. + +Reproduce the measurement from the repository root: + +```bash +python development/skeleton/benchmark_teasar_mrc.py \ + --crop-size 200 \ + --pixel-size 10 \ + --scale 0 \ + --constant 70 \ + --ball-constant 140 \ + --all-settings \ + --kimimaro \ + --threads 1 8 \ + --repeats 7 \ + --warmup 2 \ + --json /tmp/teasar_options_mrc.json +``` + +## Skeleton comparison + +The table reports raw skeletons before tick removal. A spur is a degree-3 node +with a terminal arm of at most 200 Å. A real junction has three arms longer +than 200 Å. Percentages use the number of degree-3 nodes as the denominator. + +| implementation and setting | constant | vertices | components | length | degree 1 | degree 3 | degree 4 | spurs | real junctions | +| --- | ---: | ---: | ---: | ---: | ---: | ---: | ---: | ---: | ---: | +| bioimage-cpp cube/fix | 70 Å | 11,850 | 72 | 16.723 µm | 340 | 198 | 0 | 110 (55.6%) | 22 (11.1%) | +| bioimage-cpp cube/parental | 70 Å | 11,875 | 72 | 16.727 µm | 293 | 147 | 2 | 58 (39.5%) | 22 (15.0%) | +| bioimage-cpp ball/fix | 140 Å | 11,256 | 72 | 15.857 µm | 270 | 128 | 0 | 50 (39.1%) | 26 (20.3%) | +| bioimage-cpp ball/parental | 140 Å | 11,399 | 72 | 16.039 µm | 258 | 110 | 3 | 31 (28.2%) | 25 (22.7%) | +| Kimimaro ball/fix | 140 Å | 11,328 | 71 | 16.308 µm | 283 | 141 | 0 | 67 (47.5%) | 21 (14.9%) | +| Kimimaro ball/parental | 140 Å | 11,507 | 71 | 16.505 µm | 268 | 126 | 0 | 44 (34.9%) | 24 (19.0%) | + +Ball invalidation has the largest topology effect in bioimage-cpp. Relative to +the default cube/fix mode, ball/fix reduces degree-3 nodes from 198 to 128 and +spurs from 110 to 50. The number of real junctions increases from 22 to 26. + +The fixed parental field reduces branching further. It changes the ball result +from 128 to 110 degree-3 nodes and from 50 to 31 spurs. The real-junction count +remains similar at 25. The same direction appears in Kimimaro, where the +parental mode changes 141 to 126 degree-3 nodes and 67 to 44 spurs. + +The bioimage-cpp ball modes produce fewer degree-3 nodes and spurs than their +matched Kimimaro modes. They also introduce degree-4 nodes when the parental +field is active. These metrics describe topology only; they do not establish +which skeleton is more accurate. + +## Performance comparison + +Each backend receives two warmups. The benchmark then measures seven calls in +a deterministically shuffled order. The table reports the median. MRC loading, +crop extraction, and graph-statistic calculation are outside the timed region. + +| implementation and setting | one worker | eight workers | speedup | +| --- | ---: | ---: | ---: | +| bioimage-cpp cube/fix | 0.713 s | 0.219 s | 3.26x | +| bioimage-cpp cube/parental | 0.692 s | 0.210 s | 3.30x | +| bioimage-cpp ball/fix | 0.944 s | 0.272 s | 3.47x | +| bioimage-cpp ball/parental | 0.911 s | 0.276 s | 3.30x | +| Kimimaro ball/fix | 1.833 s | 0.946 s | 1.94x | +| Kimimaro ball/parental | 1.817 s | 0.958 s | 1.90x | + +Bioimage-cpp is 1.94x faster than Kimimaro for ball/fix with one worker and +3.48x faster with eight workers. The corresponding ball/parental factors are +1.99x and 3.47x. + +The parental field reduces the one-worker bioimage-cpp time by approximately +3%. It does not provide a consistent gain with eight workers on this crop. +The saved path work is small relative to the shared distance transform and +component dispatch. + +Ball invalidation remains more expensive than cube invalidation. It is 32% +slower for fix mode with one worker and 24% slower with eight workers. This is +not a primitive-level comparison. The ball uses twice the constant radius and +produces a different set of paths. + +### Optimization stages + +The baseline is the implementation measured before this optimization pass. +The queue stage keeps only the best pending invalidation distance for each +compact foreground voxel. The final stage also shares one guarded distance +transform across ordinary binary components when its memory and volume tests +pass. + +| ball setting | workers | baseline | coalesced queue | final | final change | +| --- | ---: | ---: | ---: | ---: | ---: | +| fix | 1 | 1.985 s | 1.121 s | 0.944 s | -52.4% | +| fix | 8 | 0.518 s | 0.306 s | 0.272 s | -47.5% | +| parental | 1 | 1.945 s | 1.096 s | 0.911 s | -53.1% | +| parental | 8 | 0.516 s | 0.311 s | 0.276 s | -46.6% | + +The queue change accounts for a 39.8--43.7% reduction from the baseline. The +shared distance transform gives a further 11.2--16.8% reduction relative to +the queue stage. The measured vertices, edges, and radii did not change across +these stages. + +### Profile and memory + +A one-worker profile of ball/fix selected the shared distance transform and +estimated 53.9 MB of scratch memory. The main phase totals were: + +| phase | time | +| --- | ---: | +| shared distance-transform setup | 21.8 ms | +| shared distance transform | 151.4 ms | +| compact-domain construction | 27.0 ms | +| root Dijkstra | 283.6 ms | +| rail-path Dijkstra | 167.4 ms | +| invalidation | 273.8 ms | +| measured TEASAR total | 980.9 ms | + +The invalidation queue received 5,976,824 offers and accepted 693,630 of them. +It rejected 88.4% of offers before heap insertion. Of the accepted entries, +107,755 became stale before removal. The peak heap size was 18,059 entries. +The invalidation phase fell from approximately 1.074 s to 273.8 ms in the +profile build. The shared distance transform reduced the corresponding local +component-transform total from approximately 423 ms to 151 ms. + +Fresh worker processes measured incremental process-tree peak RSS. Each value +is the median of three processes after imports and input allocation: + +| implementation and setting | one worker | eight workers | +| --- | ---: | ---: | +| bioimage-cpp ball/fix | 52.7 MiB | 77.1 MiB | +| bioimage-cpp ball/parental | 52.7 MiB | 77.1 MiB | +| Kimimaro ball/fix | 89.8 MiB | 1,575 MiB | +| Kimimaro ball/parental | 89.8 MiB | 1,555 MiB | + +The shared transform raises the one-worker ball/fix peak from the 42.5 MiB +baseline to 52.7 MiB. It reduces the eight-worker peak from 129.8 MiB to +77.1 MiB because workers no longer retain separate component transforms. +The automatic strategy uses local transforms when the shared volume or the +estimated scratch memory exceeds its guard. + +The results cover one crop on one host. They do not predict full-volume runtime +or performance on masks with a different component-size distribution. + +## Later queue evaluation + +A follow-up tested 4-ary and radix queues for compact Dijkstra and an indexed +queue for ball invalidation. None met the retention gates, so the production +algorithm and the results in this report remain unchanged. + +The Dijkstra candidates changed a physical field by at most 1.9% and regressed +the parental field. The indexed invalidation queue removed 107,755 stale +entries on ball/fix, but one-worker MRC time changed from 974.9 ms to 982.0 ms. +Ball/parental changed from 956.6 ms to 975.8 ms. See +`PERFORMANCE_NOTES.md` for the complete queue, synthetic, and memory results. diff --git a/include/bioimage_cpp/skeleton/detail/compact_grid_dijkstra.hxx b/include/bioimage_cpp/skeleton/detail/compact_grid_dijkstra.hxx index 993b2c1..13d9268 100644 --- a/include/bioimage_cpp/skeleton/detail/compact_grid_dijkstra.hxx +++ b/include/bioimage_cpp/skeleton/detail/compact_grid_dijkstra.hxx @@ -206,14 +206,32 @@ struct CompactHeapGreater { }; struct CompactDijkstraStats { + std::size_t calls = 0; + std::size_t initialized_nodes = 0; + std::size_t settled_nodes = 0; + std::size_t neighbor_probes = 0; + std::size_t successful_relaxations = 0; + std::size_t target_marks = 0; std::size_t pushes = 0; std::size_t pops = 0; + std::size_t stale_pops = 0; std::size_t peak_heap = 0; void reset() noexcept { - pushes = 0; - pops = 0; - peak_heap = 0; + *this = {}; + } + + void merge(const CompactDijkstraStats &other) noexcept { + calls += other.calls; + initialized_nodes += other.initialized_nodes; + settled_nodes += other.settled_nodes; + neighbor_probes += other.neighbor_probes; + successful_relaxations += other.successful_relaxations; + target_marks += other.target_marks; + pushes += other.pushes; + pops += other.pops; + stale_pops += other.stale_pops; + peak_heap = std::max(peak_heap, other.peak_heap); } }; @@ -336,7 +354,8 @@ inline void compact_physical_distance_field( workspace.state.assign(n, 0); workspace.heap.clear(); if (stats != nullptr) { - stats->reset(); + ++stats->calls; + stats->initialized_nodes += n; } distances[source] = Distance{0}; workspace.state[source] = kCompactDiscovered; @@ -346,12 +365,21 @@ inline void compact_physical_distance_field( const auto entry = compact_heap_pop(workspace, stats); const auto node = entry.node; if ((workspace.state[node] & kCompactSettled) != 0) { + if (stats != nullptr) { + ++stats->stale_pops; + } continue; } workspace.state[node] |= kCompactSettled; + if (stats != nullptr) { + ++stats->settled_nodes; + } for_each_compact_neighbor( domain, node, [&](const std::uint32_t target, const double physical_length) { + if (stats != nullptr) { + ++stats->neighbor_probes; + } if ((workspace.state[target] & kCompactSettled) != 0) { return; } @@ -365,6 +393,9 @@ inline void compact_physical_distance_field( } workspace.state[target] |= kCompactDiscovered; distances[target] = candidate; + if (stats != nullptr) { + ++stats->successful_relaxations; + } compact_heap_push(workspace, {candidate, target}, stats); } ); @@ -396,7 +427,8 @@ inline void compact_node_cost_parental_field( workspace.state.assign(n, 0); workspace.heap.clear(); if (stats != nullptr) { - stats->reset(); + ++stats->calls; + stats->initialized_nodes += n; } workspace.state[source] = kCompactDiscovered; predecessors[source] = source; @@ -406,18 +438,30 @@ inline void compact_node_cost_parental_field( const auto entry = compact_heap_pop(workspace, stats); const auto node = entry.node; if ((workspace.state[node] & kCompactSettled) != 0) { + if (stats != nullptr) { + ++stats->stale_pops; + } continue; } workspace.state[node] |= kCompactSettled; + if (stats != nullptr) { + ++stats->settled_nodes; + } for_each_compact_neighbor( domain, node, [&](const std::uint32_t target, const double) { + if (stats != nullptr) { + ++stats->neighbor_probes; + } if ((workspace.state[target] & (kCompactDiscovered | kCompactSettled)) != 0) { return; } workspace.state[target] |= kCompactDiscovered; predecessors[target] = node; + if (stats != nullptr) { + ++stats->successful_relaxations; + } compact_heap_push( workspace, {static_cast(entry.distance + costs[target]), target}, @@ -459,7 +503,9 @@ inline void compact_node_cost_path( workspace.predecessors.resize(n); workspace.heap.clear(); if (stats != nullptr) { - stats->reset(); + ++stats->calls; + stats->initialized_nodes += n; + stats->target_marks += targets.size(); } for (const auto target : targets) { workspace.state[target] |= kCompactTarget; @@ -473,9 +519,15 @@ inline void compact_node_cost_path( const auto entry = compact_heap_pop(workspace, stats); const auto node = entry.node; if ((workspace.state[node] & kCompactSettled) != 0) { + if (stats != nullptr) { + ++stats->stale_pops; + } continue; } workspace.state[node] |= kCompactSettled; + if (stats != nullptr) { + ++stats->settled_nodes; + } if ((workspace.state[node] & kCompactTarget) != 0) { reached = node; break; @@ -483,6 +535,9 @@ inline void compact_node_cost_path( for_each_compact_neighbor( domain, node, [&](const std::uint32_t target, const double) { + if (stats != nullptr) { + ++stats->neighbor_probes; + } if ((workspace.state[target] & (kCompactDiscovered | kCompactSettled)) != 0) { return; @@ -492,6 +547,9 @@ inline void compact_node_cost_path( ); workspace.state[target] |= kCompactDiscovered; workspace.predecessors[target] = node; + if (stats != nullptr) { + ++stats->successful_relaxations; + } compact_heap_push(workspace, {candidate, target}, stats); } ); diff --git a/include/bioimage_cpp/skeleton/teasar.hxx b/include/bioimage_cpp/skeleton/teasar.hxx index 33b3085..adf6884 100644 --- a/include/bioimage_cpp/skeleton/teasar.hxx +++ b/include/bioimage_cpp/skeleton/teasar.hxx @@ -548,6 +548,7 @@ inline LatticeSkeletonGraph teasar_compact_impl( const TeasarOptions &options, const bool report_profile, detail::CompactBallInvalidationStats *ball_invalidation_stats = nullptr, + detail::CompactDijkstraStats *dijkstra_stats = nullptr, bioimage_cpp::detail::ActiveProfiler *component_profile = nullptr ) { if (prepared == nullptr) { @@ -722,7 +723,7 @@ inline LatticeSkeletonGraph teasar_compact_impl( } else { std::vector first_field; detail::compact_physical_distance_field( - domain, 0, dijkstra_workspace, first_field + domain, 0, dijkstra_workspace, first_field, dijkstra_stats ); Distance farthest_distance = Distance{-1}; for (std::uint32_t node = 0; node < domain.size(); ++node) { @@ -738,7 +739,7 @@ inline LatticeSkeletonGraph teasar_compact_impl( } } detail::compact_physical_distance_field( - domain, root, dijkstra_workspace, root_field + domain, root, dijkstra_workspace, root_field, dijkstra_stats ); for (const auto distance : root_field) { if (!std::isfinite(distance)) { @@ -812,7 +813,8 @@ inline LatticeSkeletonGraph teasar_compact_impl( if (!options.fix_branching) { BIOIMAGE_PROFILE_SCOPE(profile, "parental_field") detail::compact_node_cost_parental_field( - domain, root, pdrf, dijkstra_workspace, fixed_predecessors + domain, root, pdrf, dijkstra_workspace, fixed_predecessors, + dijkstra_stats ); } @@ -827,7 +829,8 @@ inline LatticeSkeletonGraph teasar_compact_impl( if (options.fix_branching) { BIOIMAGE_PROFILE_SCOPE(profile, "path_dijkstra") detail::compact_node_cost_path( - domain, target, skeleton_nodes, pdrf, dijkstra_workspace, path + domain, target, skeleton_nodes, pdrf, dijkstra_workspace, path, + dijkstra_stats ); } else { BIOIMAGE_PROFILE_SCOPE(profile, "path_from_parents") @@ -1066,12 +1069,13 @@ inline LatticeSkeletonGraph teasar_compact_prepared( detail::PreparedTeasarComponent prepared, const TeasarOptions &options, detail::CompactBallInvalidationStats *ball_invalidation_stats = nullptr, + detail::CompactDijkstraStats *dijkstra_stats = nullptr, bioimage_cpp::detail::ActiveProfiler *component_profile = nullptr ) { const ConstArrayView unused{}; return teasar_compact_impl( unused, &prepared, options, false, ball_invalidation_stats, - component_profile + dijkstra_stats, component_profile ); } @@ -1588,6 +1592,7 @@ std::vector skeletonize_components( const detail::OpenBlockFaces *open_faces = nullptr, std::vector> *precomputed_component_dbf = nullptr, std::vector *ball_invalidation_stats = nullptr, + std::vector *dijkstra_stats = nullptr, bioimage_cpp::detail::ActiveProfiler *component_profile = nullptr ) { const auto count = components.components.size(); @@ -1598,6 +1603,9 @@ std::vector skeletonize_components( if (ball_invalidation_stats != nullptr) { ball_invalidation_stats->assign(count, {}); } + if (dijkstra_stats != nullptr) { + dijkstra_stats->assign(count, {}); + } std::vector component_profiles; if (component_profile != nullptr) { component_profiles.resize(count); @@ -1656,6 +1664,8 @@ std::vector skeletonize_components( std::move(prepared), local_options, ball_invalidation_stats == nullptr ? nullptr : &ball_invalidation_stats->at(component_id), + dijkstra_stats == nullptr + ? nullptr : &dijkstra_stats->at(component_id), component_profile == nullptr ? nullptr : &component_profiles[component_id] ); @@ -1758,13 +1768,17 @@ inline SkeletonGraph teasar_with_component_edt( } std::vector results; std::vector ball_stats; + std::vector dijkstra_stats; bioimage_cpp::detail::ActiveProfiler component_profile; #ifdef BIOIMAGE_PROFILE auto *ball_stats_output = &ball_stats; + auto *dijkstra_stats_output = &dijkstra_stats; auto *component_profile_output = &component_profile; #else auto *ball_stats_output = static_cast *>(nullptr); + auto *dijkstra_stats_output = + static_cast *>(nullptr); auto *component_profile_output = static_cast(nullptr); #endif @@ -1773,7 +1787,7 @@ inline SkeletonGraph teasar_with_component_edt( results = detail_teasar::skeletonize_components( components, options, false, nullptr, nullptr, shared_edt.selected() ? &shared_edt.component_dbf : nullptr, - ball_stats_output, component_profile_output + ball_stats_output, dijkstra_stats_output, component_profile_output ); } std::vector component_ids(results.size()); @@ -1792,6 +1806,10 @@ inline SkeletonGraph teasar_with_component_edt( for (const auto &stats : ball_stats) { combined_ball_stats.merge(stats); } + detail::CompactDijkstraStats combined_dijkstra_stats; + for (const auto &stats : dijkstra_stats) { + combined_dijkstra_stats.merge(stats); + } std::fprintf( stderr, "[bioimage TEASAR diagnostics]\n" @@ -1802,7 +1820,17 @@ inline SkeletonGraph teasar_with_component_edt( " ball_pops %zu\n" " ball_stale_pops %zu\n" " ball_invalidated %zu\n" - " ball_peak_heap %zu\n", + " ball_peak_heap %zu\n" + " dijkstra_calls %zu\n" + " dijkstra_initialized %zu\n" + " dijkstra_settled %zu\n" + " dijkstra_probes %zu\n" + " dijkstra_relaxations %zu\n" + " dijkstra_target_marks %zu\n" + " dijkstra_pushes %zu\n" + " dijkstra_pops %zu\n" + " dijkstra_stale_pops %zu\n" + " dijkstra_peak_heap %zu\n", detail_teasar::shared_edt_decision_name(shared_edt.decision), shared_edt.estimated_scratch_bytes, combined_ball_stats.offers, @@ -1810,7 +1838,17 @@ inline SkeletonGraph teasar_with_component_edt( combined_ball_stats.pops, combined_ball_stats.stale_pops, combined_ball_stats.invalidated, - combined_ball_stats.peak_heap + combined_ball_stats.peak_heap, + combined_dijkstra_stats.calls, + combined_dijkstra_stats.initialized_nodes, + combined_dijkstra_stats.settled_nodes, + combined_dijkstra_stats.neighbor_probes, + combined_dijkstra_stats.successful_relaxations, + combined_dijkstra_stats.target_marks, + combined_dijkstra_stats.pushes, + combined_dijkstra_stats.pops, + combined_dijkstra_stats.stale_pops, + combined_dijkstra_stats.peak_heap ); #endif return detail_teasar::lattice_to_physical(