Break tied connected-component peaks toward the centroid

spotfinder/CMakeLists.txt

@@ -75,4 +75,6 @@ set_property(TARGET spotfinder32 PROPERTY CUDA_SEPARABLE_COMPILATION ON)
 target_compile_options(spotfinder32 PRIVATE "$<$<AND:$<CONFIG:Debug>,$<COMPILE_LANGUAGE:CUDA>>:-G>")
 target_compile_options(spotfinder32 PRIVATE "$<$<AND:$<COMPILE_LANGUAGE:CUDA>,$<OR:$<CONFIG:Debug>,$<CONFIG:RelWithDebInfo>>>:--generate-line-info>")
 
-install(TARGETS spotfinder32 RUNTIME COMPONENT Runtime)
+install(TARGETS spotfinder32 RUNTIME COMPONENT Runtime)
+
+add_subdirectory(tests)

spotfinder/connected_components/connected_components.cc

@@ -141,7 +141,30 @@ void ConnectedComponents::generate_boxes(const uint32_t width,
 
 #pragma region Reflection3D
 bool Reflection3D::is_signal_preferred(const Signal &candidate,
-                                       const Signal &current) const {
+                                       const Signal &current,
+                                       double com_x,
+                                       double com_y,
+                                       double com_z) const {
+    // Prefer the pixel nearest the centroid (see
+    // https://github.com/dials/dials/issues/3014). Squared distance
+    // avoids a per-pixel sqrt; the 0.5 offset puts the position at the
+    // pixel centre, as in center_of_mass().
+    auto distance_sq = [](const Signal &s, double cx, double cy, double cz) {
+        double dx = (static_cast<double>(s.x) + 0.5) - cx;
+        double dy = (static_cast<double>(s.y) + 0.5) - cy;
+        double dz = (static_cast<double>(s.z.value()) + 0.5) - cz;
+        return dx * dx + dy * dy + dz * dz;
+    };
+
+    double candidate_dist_sq = distance_sq(candidate, com_x, com_y, com_z);
+    double current_dist_sq = distance_sq(current, com_x, com_y, com_z);
+    if (candidate_dist_sq != current_dist_sq) {
+        return candidate_dist_sq < current_dist_sq;
+    }
+
+    // Equal distance to the centroid: fall back to a deterministic z, y, x
+    // ordering so the choice is independent of signal iteration order.
+
     // Compare z-coordinates first
     if (candidate.z.value() != current.z.value()) {
         return candidate.z.value() < current.z.value();

spotfinder/connected_components/connected_components.hpp

@@ -117,6 +117,11 @@ class Reflection3D {
         // logger.debug("Finding peak signal for reflection with {} pixels",
         //              signals_.size());
 
+        // Get the center of mass up front; used for nearest-centroid
+        // tie-breaking and for the final peak-centroid distance below.
+        auto [com_x, com_y, com_z] = center_of_mass();
+        // logger.debug("Center of mass: ({:.3f}, {:.3f}, {:.3f})", com_x, com_y, com_z);
+
         // Find the signal with the highest intensity
         const Signal *peak_signal = nullptr;
         double max_intensity = std::numeric_limits<double>::min();
@@ -151,8 +156,9 @@ class Reflection3D {
                 continue;
             }
 
-            // Deterministic tie-breaking using coordinate comparison
-            bool should_update_tie = is_signal_preferred(signal, *peak_signal);
+            // Deterministic tie-breaking: prefer the pixel nearest the centroid
+            bool should_update_tie =
+              is_signal_preferred(signal, *peak_signal, com_x, com_y, com_z);
 
             // logger.trace(
             //   "Tie at intensity {}: current ({}, {}, {}) vs peak ({}, {}, {}), "
@@ -186,10 +192,6 @@ class Reflection3D {
         //   peak_signal->linear_index,
         //   candidates_with_max_intensity);
 
-        // Get the cached or computed center of mass
-        auto [com_x, com_y, com_z] = center_of_mass();
-        // logger.debug("Center of mass: ({:.3f}, {:.3f}, {:.3f})", com_x, com_y, com_z);
-
         // Calculate the Euclidean distance
         float dx = (static_cast<float>(peak_signal->x) + 0.5f) - com_x;
         float dy = (static_cast<float>(peak_signal->y) + 0.5f) - com_y;
@@ -257,14 +259,27 @@ class Reflection3D {
     mutable std::tuple<float, float, float> com_cache_;
 
     /**
-     * @brief Determines if the first signal should be preferred over the second
-     *        in case of intensity ties using coordinate-based tie-breaking.
-     * 
+     * @brief Determines if the first signal
+     *        should be preferred over the second in case of intensity
+     *        ties.
+     *
+     * The tie is broken in favour of the pixel nearest the centroid
+     * (see https://github.com/dials/dials/issues/3014); pixels
+     * equidistant from the centroid fall back to a deterministic z, y,
+     * x ordering.
+     *
      * @param candidate The signal to compare
      * @param current The current preferred signal
+     * @param com_x The x-coordinate of the center of mass
+     * @param com_y The y-coordinate of the center of mass
+     * @param com_z The z-coordinate of the center of mass
      * @return true if candidate should be preferred over current
      */
-    bool is_signal_preferred(const Signal &candidate, const Signal &current) const;
+    bool is_signal_preferred(const Signal &candidate,
+                             const Signal &current,
+                             double com_x,
+                             double com_y,
+                             double com_z) const;
 };
 
 /**

spotfinder/tests/CMakeLists.txt

@@ -0,0 +1,23 @@
+include(GoogleTest)
+
+add_executable(test_connected_components
+    test_connected_components.cc
+    ${CMAKE_CURRENT_SOURCE_DIR}/../connected_components/connected_components.cc
+)
+target_include_directories(test_connected_components PRIVATE
+    ${CMAKE_CURRENT_SOURCE_DIR}/../connected_components
+)
+target_link_libraries(test_connected_components PRIVATE
+    ffs_common
+    dx2
+    h5read
+    Boost::graph
+    GTest::gtest_main
+    Eigen3::Eigen
+    spdlog::spdlog
+    nlohmann_json::nlohmann_json
+    CUDA::cudart
+    lodepng
+)
+
+gtest_discover_tests(test_connected_components PROPERTIES LABELS "spotfinder-tests")

spotfinder/tests/test_connected_components.cc

@@ -0,0 +1,80 @@
+/**
+ * @file test_connected_components.cc
+ * @brief Unit tests for Reflection3D peak selection (connected_components).
+ *
+ * Tests tie-breaking in peak_centroid_distance(): when several pixels
+ * share the maximum intensity, the peak is the pixel nearest the
+ * centroid rather than the first in z, y, x order (see
+ * https://github.com/dials/dials/issues/3014). The selected peak is
+ * read back through the returned peak-centroid distance.
+ */
+
+#include <gtest/gtest.h>
+
+#include <cmath>
+#include <optional>
+
+#include "connected_components.hpp"
+
+namespace {
+
+/// Build a single-frame signal (z = 0) at integer pixel (x, y).
+Signal make_signal(uint32_t x, uint32_t y, pixel_t intensity) {
+    return Signal{x, y, std::optional<int>(0), intensity, 0};
+}
+
+/// Euclidean distance from a pixel centre to the given centroid.
+float distance_to_com(uint32_t x,
+                      uint32_t y,
+                      int z,
+                      float com_x,
+                      float com_y,
+                      float com_z) {
+    float dx = (static_cast<float>(x) + 0.5f) - com_x;
+    float dy = (static_cast<float>(y) + 0.5f) - com_y;
+    float dz = (static_cast<float>(z) + 0.5f) - com_z;
+    return std::sqrt(dx * dx + dy * dy + dz * dz);
+}
+
+}  // namespace
+
+// Single brightest pixel -> the peak is that pixel and the reported
+// distance is its distance to the centroid. Covers the common untied
+// path.
+TEST(ReflectionPeakSelection, UniqueMaximum) {
+    Reflection3D reflection;
+    reflection.add_signal(make_signal(0, 0, 2));
+    reflection.add_signal(make_signal(3, 0, 8));  // unique maximum
+
+    auto [com_x, com_y, com_z] = reflection.center_of_mass();
+    float expected = distance_to_com(3, 0, 0, com_x, com_y, com_z);
+
+    EXPECT_FLOAT_EQ(reflection.peak_centroid_distance(), expected);
+}
+
+// Two pixels share the maximum value. The intensity mass (and so the
+// centroid) sits in the high-x, high-y corner, so the tie breaks to the
+// pixel nearest the centroid, not the (0, 0) pixel that wins the old z,
+// y, x ordering. See https://github.com/dials/dials/issues/3014.
+TEST(ReflectionPeakSelection, TieResolvedByCentroid) {
+    Reflection3D reflection;
+
+    // Blob of intensity in the high-x, high-y corner so the centroid sits there.
+    for (uint32_t y = 3; y <= 5; ++y) {
+        for (uint32_t x = 3; x <= 5; ++x) {
+            reflection.add_signal(make_signal(x, y, 5));
+        }
+    }
+    reflection.add_signal(make_signal(0, 0, 10));  // far, z/y/x-first maximum
+    reflection.add_signal(make_signal(5, 5, 10));  // near the centroid
+
+    auto [com_x, com_y, com_z] = reflection.center_of_mass();
+    float near_distance = distance_to_com(5, 5, 0, com_x, com_y, com_z);
+    float far_distance = distance_to_com(0, 0, 0, com_x, com_y, com_z);
+
+    // Check the near pixel is the closer of the two.
+    ASSERT_LT(near_distance, far_distance);
+
+    // The near-centroid pixel is selected, not the (0, 0) corner.
+    EXPECT_FLOAT_EQ(reflection.peak_centroid_distance(), near_distance);
+}