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

@@ -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``,

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
 

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 <typename T> class VarClusterFinder {
     T threshold_;
     std::optional<T> neighbour_threshold_; // follows threshold_ when unset
     std::optional<NDArray<T, 2>> noise_map_;
+    bool noise_map_has_negative_ = false;
     // Label image template: -1 on masked pixels, 0 elsewhere
     std::optional<NDArray<int32_t, 2>> masked_labels_;
     T seed_factor_ = 5;
@@ -166,10 +169,16 @@ template <typename T> class VarClusterFinder {
     void set_noise_map(NDView<const T, 2> 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<NDArray<T, 2>> &noise_map() const { return noise_map_; }
@@ -252,7 +261,9 @@ template <typename T> 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 <typename T> 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<T, 2> 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 <typename T> 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 <class Threshold>
     void find_clusters_impl(NDView<T, 2> img, const Threshold &thr) {
         const ssize_t n_pixels = shape_[0] * shape_[1];

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<double> img) {

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

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<double, 2> 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<double, 2> 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<double, 2> img({8, 8}, 0.0);
+    add_cluster(img, 3, 3);
+    const NDArray<double, 2> 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<double, 2> 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<double, 2> 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",