diff --git a/spotfinder/CMakeLists.txt b/spotfinder/CMakeLists.txt index 4a3eb9c4..2f9c4ac2 100644 --- a/spotfinder/CMakeLists.txt +++ b/spotfinder/CMakeLists.txt @@ -75,4 +75,6 @@ set_property(TARGET spotfinder32 PROPERTY CUDA_SEPARABLE_COMPILATION ON) target_compile_options(spotfinder32 PRIVATE "$<$,$>:-G>") target_compile_options(spotfinder32 PRIVATE "$<$,$,$>>:--generate-line-info>") -install(TARGETS spotfinder32 RUNTIME COMPONENT Runtime) \ No newline at end of file +install(TARGETS spotfinder32 RUNTIME COMPONENT Runtime) + +add_subdirectory(tests) \ No newline at end of file diff --git a/spotfinder/connected_components/connected_components.cc b/spotfinder/connected_components/connected_components.cc index e8354d03..c859de31 100644 --- a/spotfinder/connected_components/connected_components.cc +++ b/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 ¤t) const { + const Signal ¤t, + 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(s.x) + 0.5) - cx; + double dy = (static_cast(s.y) + 0.5) - cy; + double dz = (static_cast(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(); diff --git a/spotfinder/connected_components/connected_components.hpp b/spotfinder/connected_components/connected_components.hpp index ca896a15..cbd5b9ec 100644 --- a/spotfinder/connected_components/connected_components.hpp +++ b/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::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(peak_signal->x) + 0.5f) - com_x; float dy = (static_cast(peak_signal->y) + 0.5f) - com_y; @@ -257,14 +259,27 @@ class Reflection3D { mutable std::tuple 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 ¤t) const; + bool is_signal_preferred(const Signal &candidate, + const Signal ¤t, + double com_x, + double com_y, + double com_z) const; }; /** diff --git a/spotfinder/tests/CMakeLists.txt b/spotfinder/tests/CMakeLists.txt new file mode 100644 index 00000000..8dc1a461 --- /dev/null +++ b/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") diff --git a/spotfinder/tests/test_connected_components.cc b/spotfinder/tests/test_connected_components.cc new file mode 100644 index 00000000..2ac0f671 --- /dev/null +++ b/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 + +#include +#include + +#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(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(x) + 0.5f) - com_x; + float dy = (static_cast(y) + 0.5f) - com_y; + float dz = (static_cast(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); +}