From 48184ff71fabe600cf7c2d46382aa0eeaa3e6b8d Mon Sep 17 00:00:00 2001 From: kferjaoui Date: Tue, 3 Nov 2026 17:26:31 +0100 Subject: [PATCH] Require 0 <= neighbour <= seed threshold when erasing clusters With erase_clusters on, visit() zeroes every neighbour that fails the neighbour test. If the neighbour threshold was above the seed threshold, that could zero a pixel that would later have seeded its own cluster. With a negative threshold the zero written back could itself pass a test and create seeds or cluster pixels. --- RELEASE.md | 5 +- .../src/python/cluster/pyVarClusterFinder.rst | 7 +- include/aare/VarClusterFinder.hpp | 51 ++++++++++++-- python/src/var_cluster.hpp | 5 +- python/tests/test_VarClusterFinder.py | 16 +++++ src/VarClusterFinder.test.cpp | 67 +++++++++++++++++++ 6 files changed, 142 insertions(+), 9 deletions(-) diff --git a/RELEASE.md b/RELEASE.md index f7018a18..839fc092 100644 --- a/RELEASE.md +++ b/RELEASE.md @@ -40,7 +40,10 @@ modified unless erasing is enabled (``set_erase_clusters`` in C++, the ``erase_clusters`` property in Python, replacing ``set_empty_surroundingPixels``), which now zeroes both the cluster - pixels and their sub-threshold neighbours. + pixels and their sub-threshold neighbours. Erasing requires + 0 <= neighbour threshold <= seed threshold (with a noise map, ordered + non-negative factors and no negative noise) so that it never changes + the hits; ``find_clusters`` raises an error otherwise. Hits are stored in an offset layout instead of one 624 byte record per cluster: ``steal_hits`` returns a ``Hits`` structure (``VarClusterHits`` in Python) with one 32 byte record per cluster (``size``, ``row``, diff --git a/docs/src/python/cluster/pyVarClusterFinder.rst b/docs/src/python/cluster/pyVarClusterFinder.rst index 2d4ce4c7..2ed74452 100644 --- a/docs/src/python/cluster/pyVarClusterFinder.rst +++ b/docs/src/python/cluster/pyVarClusterFinder.rst @@ -95,8 +95,11 @@ the constructor and leaves it untouched. With ``erase_clusters = True`` the clusters are erased from the frame as the fill visits them: every cluster pixel is set to zero, and so are the sub-threshold neighbours around it, because charge that leaked below threshold would otherwise bias a pedestal -updated from the erased frame. A zeroed pixel compares as zero for the rest -of the frame, which only matters for negative thresholds. +updated from the erased frame. Erasing requires +0 <= ``neighbour_threshold`` <= ``threshold`` (with a noise map, +0 <= ``neighbour_factor`` <= ``seed_factor`` and no negative noise), so that +it never changes the hits; ``find_clusters`` raises ``RuntimeError`` +otherwise. .. py:currentmodule:: aare diff --git a/include/aare/VarClusterFinder.hpp b/include/aare/VarClusterFinder.hpp index 62da9131..643c5ba3 100644 --- a/include/aare/VarClusterFinder.hpp +++ b/include/aare/VarClusterFinder.hpp @@ -46,8 +46,10 @@ namespace aare { * in which case every cluster pixel, and also every sub-threshold neighbour * the fill examined around it, is set to zero as it is visited. The ring is * included so that a pedestal updated from the erased frame sees no charge - * that leaked out of the cluster. A zeroed pixel compares as zero for the - * rest of the frame, which only matters for negative thresholds. + * that leaked out of the cluster. Erasing requires 0 <= neighbour threshold + * <= seed threshold at every pixel (with a noise map: 0 <= neighbour factor + * <= seed factor and no negative noise), so that it never changes the hits; + * find_clusters() throws otherwise. * * @tparam T Pixel type of the image and the noise map. */ @@ -118,6 +120,7 @@ template class VarClusterFinder { T threshold_; std::optional neighbour_threshold_; // follows threshold_ when unset std::optional> noise_map_; + bool noise_map_has_negative_ = false; // Label image template: -1 on masked pixels, 0 elsewhere std::optional> masked_labels_; T seed_factor_ = 5; @@ -166,10 +169,16 @@ template class VarClusterFinder { void set_noise_map(NDView noise_map) { check_shape(noise_map.shape(), "noise map"); noise_map_.emplace(noise_map); + noise_map_has_negative_ = + std::any_of(noise_map_->begin(), noise_map_->end(), + [](T v) { return v < T{}; }); } /** @brief Drop the noise map and use the fixed threshold again. */ - void clear_noise_map() { noise_map_.reset(); } + void clear_noise_map() { + noise_map_.reset(); + noise_map_has_negative_ = false; + } /** @brief The noise map, or nothing while the fixed threshold is used. */ const std::optional> &noise_map() const { return noise_map_; } @@ -252,7 +261,9 @@ template class VarClusterFinder { * @brief Erase the clusters from the input image while processing a * frame (default off): every cluster pixel and the sub-threshold * neighbours around it are set to zero, so that a pedestal updated from - * the frame afterwards is not biased by the clusters. + * the frame afterwards is not biased by the clusters. Requires + * 0 <= neighbour threshold <= seed threshold; find_clusters() throws + * otherwise. */ void set_erase_clusters(bool erase) { erase_clusters_ = erase; } /** @brief Whether clusters are erased from the input image. */ @@ -261,10 +272,14 @@ template class VarClusterFinder { /** * @brief Find the clusters in one frame and append them to the hits. * @param img Image of the shape given to the constructor. - * @throws std::runtime_error if the shape does not match. + * @throws std::runtime_error if the shape does not match, or if erasing + * is enabled and the thresholds do not satisfy + * 0 <= neighbour threshold <= seed threshold (see set_erase_clusters()). */ void find_clusters(NDView img) { check_shape(img.shape(), "image"); + if (erase_clusters_) + check_erase_thresholds(); if (noise_map_) find_clusters_impl(img, NoiseThreshold{noise_map_->data(), seed_factor_, @@ -298,6 +313,32 @@ template class VarClusterFinder { what, shape[0], shape[1], shape_[0], shape_[1])); } + // Erasing must not change the hits: a zeroed pixel has to fail both + // tests, and a pixel rejected as a neighbour must not pass as a seed. + // Both hold when 0 <= neighbour(i) <= seed(i) for every pixel. + void check_erase_thresholds() const { + if (noise_map_) { + if (!(0 <= neighbour_factor_ && neighbour_factor_ <= seed_factor_)) + throw std::runtime_error( + LOCATION + + fmt::format("erase_clusters requires 0 <= neighbour " + "factor ({}) <= seed factor ({})", + neighbour_factor_, seed_factor_)); + if (noise_map_has_negative_) + throw std::runtime_error(LOCATION + + "erase_clusters requires a noise map " + "without negative values"); + } else { + const T neighbour = neighbour_threshold(); + if (!(T{} <= neighbour && neighbour <= threshold_)) + throw std::runtime_error( + LOCATION + + fmt::format("erase_clusters requires 0 <= neighbour " + "threshold ({}) <= threshold ({})", + neighbour, threshold_)); + } + } + template void find_clusters_impl(NDView img, const Threshold &thr) { const ssize_t n_pixels = shape_[0] * shape_[1]; diff --git a/python/src/var_cluster.hpp b/python/src/var_cluster.hpp index 9aade712..505c89a1 100644 --- a/python/src/var_cluster.hpp +++ b/python/src/var_cluster.hpp @@ -189,7 +189,10 @@ per pixel multiples of the noise, ``seed_factor`` and ``neighbour_factor``.)") "Erase the clusters from the input image while processing a " "frame (default False): every cluster pixel and the sub-threshold " "neighbours around it are set to zero, so that a pedestal updated " - "from the frame afterwards is not biased by the clusters.") + "from the frame afterwards is not biased by the clusters. " + "Requires 0 <= neighbour threshold <= threshold (with a noise " + "map: 0 <= neighbour_factor <= seed_factor and no negative " + "noise); find_clusters raises RuntimeError otherwise.") .def( "find_clusters", [](Finder &self, py::array_t img) { diff --git a/python/tests/test_VarClusterFinder.py b/python/tests/test_VarClusterFinder.py index fd475f37..0968ce69 100644 --- a/python/tests/test_VarClusterFinder.py +++ b/python/tests/test_VarClusterFinder.py @@ -134,6 +134,22 @@ def test_erase_clusters_zeroes_cluster_and_ring(): assert np.array_equal(img, expected) +def test_erase_clusters_requires_ordered_non_negative_thresholds(): + img = frame_with_cluster() + cf = VarClusterFinder((8, 8), 10.0) + cf.erase_clusters = True + cf.neighbour_threshold = 30.0 + with pytest.raises(RuntimeError, match="erase_clusters requires"): + cf.find_clusters(img) + assert img[3, 4] == 100.0 # rejected before touching anything + cf.neighbour_threshold = None + noise = np.ones((8, 8)) + noise[0, 0] = -1.0 + cf.noise_map = noise + with pytest.raises(RuntimeError, match="negative"): + cf.find_clusters(img) + + def test_labeled_image(): img = frame_with_cluster() img[7, 0] = 50.0 diff --git a/src/VarClusterFinder.test.cpp b/src/VarClusterFinder.test.cpp index 06c9b75b..d1fa462b 100644 --- a/src/VarClusterFinder.test.cpp +++ b/src/VarClusterFinder.test.cpp @@ -395,6 +395,73 @@ TEST_CASE("The input is left untouched unless erase_clusters is enabled", } CHECK(copy(3, 5) == 0.0); } + SECTION("Zeroing does not change the hits with a lower neighbour " + "threshold") { + NDArray two({8, 12}, 1.0); + add_cluster(two, 3, 3); + add_cluster(two, 3, 7); + two(1, 3) = 7.0; // joins the first cluster as a neighbour only + NDArray copy(two); + + Finder plain({8, 12}, 10.0); + plain.set_neighbour_threshold(5.0); + plain.find_clusters(two.view()); + Finder zeroing({8, 12}, 10.0); + zeroing.set_neighbour_threshold(5.0); + zeroing.set_erase_clusters(true); + zeroing.find_clusters(copy.view()); + + const auto a = plain.steal_hits(); + const auto b = zeroing.steal_hits(); + REQUIRE(a.size() == 2); + REQUIRE(b.size() == 2); + CHECK(a.hits[0].size == 10); + for (size_t i = 0; i < a.size(); ++i) { + CHECK(a.hits[i].size == b.hits[i].size); + CHECK(a.hits[i].energy == b.hits[i].energy); + CHECK(pixels_of(a, i) == pixels_of(b, i)); + } + } +} + +TEST_CASE("Erasing requires 0 <= neighbour threshold <= seed threshold", + "[varcluster]") { + NDArray img({8, 8}, 0.0); + add_cluster(img, 3, 3); + const NDArray original(img); + Finder finder({8, 8}, 10.0); + finder.set_erase_clusters(true); + + SECTION("Neighbour threshold above the seed threshold") { + finder.set_neighbour_threshold(30.0); + CHECK_THROWS_AS(finder.find_clusters(img.view()), std::runtime_error); + CHECK(img == original); // rejected before touching anything + CHECK(finder.total_clusters() == 0); + finder.set_erase_clusters(false); + CHECK_NOTHROW(finder.find_clusters(img.view())); + } + SECTION("Negative seed threshold") { + finder.set_threshold(-5.0); + CHECK_THROWS_AS(finder.find_clusters(img.view()), std::runtime_error); + } + SECTION("Negative neighbour threshold") { + finder.set_neighbour_threshold(-1.0); + CHECK_THROWS_AS(finder.find_clusters(img.view()), std::runtime_error); + } + SECTION("Neighbour factor above the seed factor") { + NDArray noise({8, 8}, 1.0); + finder.set_noise_map(noise.view()); + finder.set_neighbour_factor(6.0); + CHECK_THROWS_AS(finder.find_clusters(img.view()), std::runtime_error); + } + SECTION("Negative noise") { + NDArray noise({8, 8}, 1.0); + noise(7, 7) = -1.0; + finder.set_noise_map(noise.view()); + CHECK_THROWS_AS(finder.find_clusters(img.view()), std::runtime_error); + finder.clear_noise_map(); + CHECK_NOTHROW(finder.find_clusters(img.view())); + } } TEST_CASE("Masked pixels never seed and are never added to a cluster",