From 8f1b0b2281438d2181e50e93c43d38efb1d45b70 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 2 Aug 2026 19:03:35 +0200 Subject: [PATCH] spot_finding: accumulate spot centroids in integers The photon-weighted position sums were floats, so the centroid's last bit depended on the build rather than on the data: gcc contracts the multiply-add in AddPixel into an FMA under -march=x86-64-v3 and cannot at the baseline, and MSVC does not contract at all under /fp:precise. The GPU extractor had to match with __fmaf_rn, and the parity test still needed a two-ulp slack for hosts that do not fuse. Column, line and the per-pixel count are all integral, so the sums are exact in int64 and both implementations reach the same bits with nothing to match. The parity test now demands exact equality unconditionally and gets it, including on a baseline build. ConvertToImageCoordinates keeps the sums integral too: the raw -> image map is a signed axis swap plus an integer translation, so it is applied to the sums instead of to the centroid. Drops the SpotToSave constructor, which had no callers and could not have been converted without quantising the stored centroid. Co-Authored-By: Claude Opus 5 (1M context) --- common/DiffractionSpot.cpp | 49 +++++++++++-------- common/DiffractionSpot.h | 18 ++++--- .../spot_finding/SpotExtractorGPU.cu | 16 +++--- .../spot_finding/SpotExtractorGPU.h | 10 ++-- tests/SpotExtractorGPUParityTest.cpp | 30 ++---------- 5 files changed, 56 insertions(+), 67 deletions(-) diff --git a/common/DiffractionSpot.cpp b/common/DiffractionSpot.cpp index bb0c473f..ba58bbe2 100644 --- a/common/DiffractionSpot.cpp +++ b/common/DiffractionSpot.cpp @@ -7,8 +7,8 @@ DiffractionSpot::DiffractionSpot(uint32_t col, uint32_t line, int64_t in_photons) { if (in_photons < 0) in_photons = 0; - x = col * static_cast(in_photons); - y = line * static_cast(in_photons); + x = static_cast(col) * in_photons; + y = static_cast(line) * in_photons; pixel_count = 1; photons = in_photons; max_photons = in_photons; @@ -23,18 +23,10 @@ DiffractionSpot& DiffractionSpot::operator+=(const DiffractionSpot &other) { return *this; } -DiffractionSpot::DiffractionSpot(float x_sum, float y_sum, int64_t in_pixel_count, +DiffractionSpot::DiffractionSpot(int64_t x_sum, int64_t y_sum, int64_t in_pixel_count, int64_t in_photons, int64_t in_max_photons) : x(x_sum), y(y_sum), pixel_count(in_pixel_count), photons(in_photons), max_photons(in_max_photons) {} -DiffractionSpot::DiffractionSpot(const SpotToSave &save) { - x = save.x * static_cast(save.intensity); - y = save.y * static_cast(save.intensity); - pixel_count = 1; - photons = std::lround(save.intensity); - max_photons = save.maxc; -} - int64_t DiffractionSpot::Count() const { return photons; } @@ -46,7 +38,10 @@ int64_t DiffractionSpot::MaxCount() const { Coord DiffractionSpot::RawCoord() const { if (photons == 0) return {0, 0, 0}; - return {x / (float)photons, y / (float)photons, 0}; + // In double: the sums run past what a float mantissa holds on a bright spot, and the centroid is + // wanted to better than the last pixel bit. + return {static_cast(static_cast(x) / static_cast(photons)), + static_cast(static_cast(y) / static_cast(photons)), 0}; } int64_t DiffractionSpot::PixelCount() const { @@ -54,24 +49,38 @@ int64_t DiffractionSpot::PixelCount() const { } void DiffractionSpot::AddPixel(uint32_t col, uint32_t line, int64_t photons) { - this->x += col * (float) photons; - this->y += line * (float) photons; + this->x += static_cast(col) * photons; + this->y += static_cast(line) * photons; this->photons += photons; this->max_photons = std::max(this->max_photons, photons); this->pixel_count += 1; } void DiffractionSpot::ConvertToImageCoordinates(const DiffractionExperiment &experiment, uint16_t module_number) { - auto c_out = RawToConvertedCoordinate(experiment, module_number, RawCoord()); - this->x = c_out.x * (float) photons; - this->y = c_out.y * (float) photons; + // The raw -> image map is a signed axis swap (module axes are +/-1 along X or Y, never rotated - + // see DetectorGeometryModular::GetDirection) plus an integer translation: the module's origin in + // the assembled image, and the two-pixel multipixel gaps, which depend on where the centroid + // falls. Applied to the SUMS rather than to the centroid, it therefore leaves them exact + // integers; converting the centroid and multiplying it back by the photon count would put the + // rounding this class avoids straight back in. + const Coord centroid = RawCoord(); + const Coord fast = experiment.GetModuleFastDirection(module_number); + const Coord slow = experiment.GetModuleSlowDirection(module_number); + // Everything the map adds on top of the axis swap - origin plus gap corrections, an integer. + const Coord shift = RawToConvertedCoordinate(experiment, module_number, centroid) + - fast * centroid.x - slow * centroid.y; + + const int64_t x_raw = x, y_raw = y; + x = std::lround(shift.x) * photons + std::lround(fast.x) * x_raw + std::lround(slow.x) * y_raw; + y = std::lround(shift.y) * photons + std::lround(fast.y) * x_raw + std::lround(slow.y) * y_raw; } std::optional DiffractionSpot::Export(const DiffractionGeometry &geometry, int64_t image_num) const { if (photons == 0) return std::nullopt; - auto d = geometry.PxlToRes(x / (float) photons, y / (float) photons); + const Coord centroid = RawCoord(); + auto d = geometry.PxlToRes(centroid.x, centroid.y); float phi = 0.0f; if (geometry.GetRotation()) { @@ -81,8 +90,8 @@ std::optional DiffractionSpot::Export(const DiffractionGeometry &geo } return SpotToSave{ - .x = x / static_cast(photons), - .y = y / static_cast(photons), + .x = centroid.x, + .y = centroid.y, .phi = phi, .intensity = static_cast(photons), .maxc = max_photons, diff --git a/common/DiffractionSpot.h b/common/DiffractionSpot.h index e98fc34b..482f988c 100644 --- a/common/DiffractionSpot.h +++ b/common/DiffractionSpot.h @@ -10,19 +10,23 @@ // Definition of Bragg spot class DiffractionSpot { - float x = 0; - float y = 0; + // Photon-weighted position sums, sum(col * photons) and sum(line * photons). Integers, not a + // centroid: column, line and the per-pixel count are all integral, so the sums are exact and + // every implementation that builds them - host, GPU, whatever the build lets the compiler + // contract - arrives at the same bits. In float they did not: gcc fuses the multiply-add under + // -march=x86-64-v3 and not at the baseline, which left the last bit of the centroid a property + // of the build flags rather than of the data. + int64_t x = 0; + int64_t y = 0; int64_t pixel_count = 0; int64_t photons = 0; // total photon count int64_t max_photons = INT64_MIN; // maximum number of counts per pixel in the spot public: DiffractionSpot() = default; DiffractionSpot(uint32_t col, uint32_t line, int64_t photons); - DiffractionSpot(const SpotToSave &save); - // From already-summed quantities. x_sum/y_sum are sum(col * photons) / sum(line * photons), i.e. - // the members as AddPixel leaves them, NOT a centroid. Used by the GPU spot extractor, which - // computes those sums on the device. - DiffractionSpot(float x_sum, float y_sum, int64_t pixel_count, int64_t photons, int64_t max_photons); + // From already-summed quantities, i.e. the members as AddPixel leaves them, NOT a centroid. + // Used by the GPU spot extractor, which computes those sums on the device. + DiffractionSpot(int64_t x_sum, int64_t y_sum, int64_t pixel_count, int64_t photons, int64_t max_photons); DiffractionSpot& operator+=(const DiffractionSpot& spot); int64_t PixelCount() const; int64_t Count() const; diff --git a/image_analysis/spot_finding/SpotExtractorGPU.cu b/image_analysis/spot_finding/SpotExtractorGPU.cu index 09a1ac7f..5c8d7c99 100644 --- a/image_analysis/spot_finding/SpotExtractorGPU.cu +++ b/image_analysis/spot_finding/SpotExtractorGPU.cu @@ -212,7 +212,7 @@ __global__ void finish_components(const uint32_t *__restrict__ index, const int3 for (int i = t; i < n; i += nthreads) atomicAdd(&count[label[root[i]]], 1); __syncthreads(); - // 3) sums, one thread per component, walking its members in ascending list order so the float + // 3) sums, one thread per component, walking its members in ascending list order so the // accumulation matches DiffractionSpot::AddPixel term for term. A component bigger than // max-pix is thrown away below, so it is not summed - which is also what keeps a whole lit // module or diffraction ring from turning into one thread walking tens of thousands of @@ -223,20 +223,16 @@ __global__ void finish_components(const uint32_t *__restrict__ index, const int3 const int want = count[l]; scratch[l].pixel_count = want; if (want > max_pix) continue; - float x = 0.0f, y = 0.0f; + long long x = 0, y = 0; long long photons = 0, max_photons = LLONG_MIN; int found = 0; for (int j = i; j < n && found < want; j++) { if (root[j] != static_cast(i)) continue; const long long counts = value[j]; - // Spelled out rather than left as "x += col * counts", because the rounding has to match - // DiffractionSpot::AddPixel term for term and the two compilers do not contract alike: - // gcc fuses AddPixel into a vfmadd whenever the build enables FMA (the CI flags do), so - // __fmaf_rn is the counterpart there. Where the host cannot fuse (a baseline -march, or - // MSVC, which does not contract by default) this leaves the last bit of the centroid - // differing by an ulp - see SpotExtractorGPUParityTest.cpp. - x = __fmaf_rn(static_cast(index[j] % width), static_cast(counts), x); - y = __fmaf_rn(static_cast(index[j] / width), static_cast(counts), y); + // Integers, exactly as DiffractionSpot::AddPixel does them, so host and device agree by + // construction - no rounding mode to match and nothing for either compiler to contract. + x += static_cast(index[j] % width) * counts; + y += static_cast(index[j] / width) * counts; photons += counts; max_photons = max(max_photons, counts); found++; diff --git a/image_analysis/spot_finding/SpotExtractorGPU.h b/image_analysis/spot_finding/SpotExtractorGPU.h index b3eafa2e..e00a86cc 100644 --- a/image_analysis/spot_finding/SpotExtractorGPU.h +++ b/image_analysis/spot_finding/SpotExtractorGPU.h @@ -21,9 +21,9 @@ // * both make a component's root its lowest list index, so both find the same roots; // * labels are handed out by a prefix sum over the roots in ascending order, which is the order the // host's second scan hands them out in, so the SPOT ORDER is identical; -// * the centroid sums are accumulated per component in ascending list order, in float, term for -// term as DiffractionSpot::AddPixel does them, with the rounding spelled out (see the comment at -// the sum itself - the two compilers do not contract a multiply-add alike). +// * the centroid sums are accumulated per component in ascending list order, in integers, term for +// term as DiffractionSpot::AddPixel does them, so there is no rounding for the two compilers to +// disagree about. // tests/SpotExtractorGPUParityTest.cpp holds the two to each other on realistic, occupancy-swept and // pathological frames, and checks that repeating a frame gives byte-identical output. @@ -38,8 +38,8 @@ // Per-component sums, in exactly the form DiffractionSpot holds them: x and y are sum(col*photons) // and sum(line*photons), not a centroid. struct SpotExtractorGPUSpot { - float x; - float y; + int64_t x; + int64_t y; int64_t photons; int64_t max_photons; int32_t pixel_count; diff --git a/tests/SpotExtractorGPUParityTest.cpp b/tests/SpotExtractorGPUParityTest.cpp index 622d30e4..b5a86d9b 100644 --- a/tests/SpotExtractorGPUParityTest.cpp +++ b/tests/SpotExtractorGPUParityTest.cpp @@ -195,26 +195,11 @@ public: } }; -// How many representable floats apart two values are. -int64_t UlpDistance(float a, float b) { - int32_t ia, ib; - memcpy(&ia, &a, sizeof(ia)); - memcpy(&ib, &b, sizeof(ib)); - if (ia < 0) ia = INT32_MIN - ia; // map to a monotone ordering across the sign - if (ib < 0) ib = INT32_MIN - ib; - return std::abs(static_cast(ia) - static_cast(ib)); -} - -// Everything that decides which spots exist and what they weigh is compared EXACTLY: the number of -// spots, their order, and each one's pixel count, photon sum and maximum. Those are integers, and a -// difference in any of them is a difference in the partition. -// -// The centroid is a float sum, and its last bit is a property of the BUILD rather than of either -// implementation: gcc contracts DiffractionSpot::AddPixel into an FMA whenever the flags allow it -// (the CI -march=x86-64-v3 does), while a baseline -march, or MSVC with its default /fp:precise, -// cannot. The extractor uses __fmaf_rn, so it is bit-exact against a host that fuses and one ulp -// off one that does not. An ulp bound catches a real divergence - which moves a centroid by pixels, -// not by 1e-4 of one - while staying true whichever way the host was built. +// Everything is compared EXACTLY: the number of spots, their order, and each one's pixel count, +// photon sum, maximum and centroid. All of them are built out of integer sums on both sides, so +// there is nothing here that a build flag can move - which is the point of accumulating in integers +// rather than in float, where gcc contracted the multiply-add under -march=x86-64-v3 and not at the +// baseline and left the last bit of the centroid a property of how the host was built. void RequireIdentical(const std::string &what, const std::vector &cpu, const std::vector &gpu) { @@ -228,13 +213,8 @@ void RequireIdentical(const std::string &what, // RawCoord divides the sums by the photon count, so comparing it compares the sums; a spot // whose photons sum to zero reports (0,0) on both sides by the same branch. const Coord c = cpu[i].RawCoord(), g = gpu[i].RawCoord(); -#ifdef __FMA__ REQUIRE(memcmp(&g.x, &c.x, sizeof(float)) == 0); REQUIRE(memcmp(&g.y, &c.y, sizeof(float)) == 0); -#else - REQUIRE(UlpDistance(g.x, c.x) <= 2); - REQUIRE(UlpDistance(g.y, c.y) <= 2); -#endif } }