Adaptive spot finding: the ring threshold replaces the floor, not the local test

The self-calibrating finder was meant to replace the classic finder's FIXED PHOTON FLOOR with a per-resolution-ring threshold read off the image's own noise. As written it replaced the local-box SNR test as well, and that is the defect: a whole-ring threshold is an ABSOLUTE contour with no feedback from a pixel's own surroundings, so the area a spot puts above it grows as sigma^2*ln(peak/threshold) and never saturates. Measured on a strongly diffracting rotation set, the detected footprint grows by +8.05 pixels per e-fold of peak, so the brightest reflections came out as 100-500 pixel blobs and were then discarded for exceeding the size bound - every one of the ten strongest on an image. Intersecting with the local box gives -0.24 pixels per e-fold, the classic finder's own number to two decimals. WHY the local box is the right partner, rather than merely the incumbent: it is a prominence rule whose reference level is a 961-pixel mean. A spot inflates the box's own variance and the peak divides out of the acceptance test, so it cuts at a fixed FRACTION of the spot's own height. Referring that level to fewer pixels makes it inherit their shot noise - at FIXED footprint, estimating the level from 961 pixels, from 25, and from the single maximum gives centroid residuals of 0.524, 0.539 and 0.656 - so flat growth and a stable centroid turn out to be two ends of one dial. A contour on the bare maximum has the flattest growth of anything tried (+0.1) and merges worst. The two arms bind in different regimes, which is why intersecting beats choosing: on serial stills the ring threshold is 0.6x the classic floor, on this rotation sweep 2.3-6.0x. Stills are a strict no-op - 175 components against 175, identical per frame - so the +40% in stills indexing that the adaptive threshold was introduced for is untouched. What it buys, stated as one fact rather than two. Across five geometry pins spanning 1.1 mm it indexes the most frames of any arm tried, 0.831 against 0.803, and integrates 3.04 to 5.76% more observations - but those are the SAME number: regressing observation count on indexing rate over four arms leaves residuals of +/-0.7 percentage points against swings of -7 to +4.5%, so the extra observations ARE the extra indexed frames, not better data per frame. CC1/2, the only statistic here carrying per-observation quality, is +0.66 at one pin and -0.06 at the other: not harmed, not improved. <I/sigma>, ISa and R_meas cannot arbitrate on this data - across those pins each crosses zero as a monotone function of the pin. WHY an absolute contour indexes fewer frames, when its spot list is equal or better on every axis measured - recall, top-1000 recall, centroid, ice fraction, component count - is the interesting part, and it is not a detection effect at all: ITS OWN SIZE BOUND DELETES THE BRIGHTEST REFLECTIONS ON THE FRAME. A component is discarded because it grew past 200 px, and it grew past 200 px because it was bright, so the deletions are drawn from the head of the indexing budget rather than uniformly from it: they are 11x enriched in the top 250 of the thousand spots handed to the indexer, and the bound's own real deletions sit at MEDIAN RANK 12. Turning the bound off recovers 66% and 50% of the deficit at the two pins, against a bar registered at 33% before the run. Three of us dismissed this for most of a day on the grounds that the gates delete only ~4% of what is detected. That arithmetic was right and the denominator was wrong - a rate is not an impact when the thing being lost is selected for the property that makes it matter. Reworking the bound instead was measured and rejected: it recovers half the deficit, and it cannot be done without re-admitting what the bound is for - 68 components past 200 px, of which 8 are real and 60 are junk, where the intersect gets the 8 without the 60. The residual once the bound is off, +1.08%/+1.70%, is the contour itself. Component merging is ruled out separately: geometrically impossible here, 33.9 px minimum reflection separation against components spanning 10 px. So is a ranking effect - the intersect's lead runs +0.06% at --max-spots 250, +3.46% at 1000 and +14.26% at 2000, which is backwards for a selection artefact. Costs 0.48 ms per image in the finder, and 0.044 px of bright-spot centroid precision - measured convention-free, by fitting a line to a reflection's own centroid across five frames, after an XDS-referenced figure proved to be four fifths aperture convention. It also makes the compactness gate above it safe. On the absolute contour that gate is net damage, deleting 37 genuine reflections per ten frames; once the footprint stops growing nothing reaches its threshold at all. Also fixes a real but unexercised defect in PoissonThreshold, where the exact tail handed over to a normal approximation with a step. It changes nothing here: the clipped ring sigma is over-dispersed 1.2-4.9x against sqrt(mu) because it still contains diffraction, so the Gaussian arm wins every ring above mu=50 and none of the 522 thresholds move. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01FBumeJVx4oeXxiBRpkrE5H
2026-08-28 15:30:39 +02:00
parent 844b461232
commit 9f2696c1ac
10 changed files with 215 additions and 61 deletions
@@ -1,8 +1,10 @@
 # Changelog
 ## 1.0.0
 ### 1.0.0-rc.166
+* Self-calibrating spot detection intersects its per-resolution-ring threshold with the local signal-to-noise test again, instead of replacing it. The ring threshold takes the place of the fixed photon floor and nothing else; standing alone it followed a bright reflection's skirt outwards, so on a strongly diffracting crystal the brightest reflections were detected as 100-500 pixel blobs and then discarded for being too large.
+* The self-calibrating threshold no longer steps where it switches from the exact Poisson tail to a normal approximation.
 * A spot may hold up to 200 connected pixels rather than 50. Under the self-calibrating threshold a reflection's footprint grows with its brightness, so the old bound acted as an intensity ceiling and discarded the strongest reflections of a strongly diffracting crystal - on one such set, every one of the ten brightest on an image.
-* A spot larger than 50 pixels must also be compact, filling a fifth of the square its bounding box fits inside. An ice arc, a cosmic-ray track or a lit detector row is refused however many pixels it has, which is what the size bound used to do and now does without capping brightness.
+* A spot larger than 50 pixels must also be compact, filling a fifth of the square its bounding box fits inside, so that an ice arc, a cosmic-ray track or a lit detector row is refused however many pixels it has. It bounds the shape of what the raised size limit admits; on data without that population it changes nothing.
 * Eight further hexagonal-ice bands between 1.472 and 1.170 A are recognised. The measured list they extend stops at 1.522 A because its source went no further, not because ice does, and on a detector reaching past 1.5 A the unlisted rings were left in the data - on a strongly diffracting set they were 44% of every image's spots. Ice handling still only applies where a run trips the ice gate.
 * The rugnux results report carries `JFJOCH_DATASET_SETTINGS=`, the geometry the run integrated at written as one line of JSON in the form `jfjoch_broker` takes it, so a refined beam centre and distance can be carried back to the instrument.
 * `/entry/MX/strongPixels` and the strong-pixel plot are filled on the CPU/GPU analysis path, not only behind the FPGA, so an image the spot finder gave up on can be told from one that did not diffract.
@@ -9,8 +9,8 @@
 #include "AdaptiveThreshold.h"

 AdaptiveSpotFinderCPU::AdaptiveSpotFinderCPU(const AzimuthalIntegrationMapping &in_mapping)
-    : ImageSpotFinder(static_cast<int32_t>(in_mapping.GetWidth()),
-                      static_cast<int32_t>(in_mapping.GetHeight())),
+    : ImageSpotFinderCPU(static_cast<int32_t>(in_mapping.GetWidth()),
+                         static_cast<int32_t>(in_mapping.GetHeight())),
      mapping(in_mapping) {
    const size_t nbins = mapping.GetBinNumber();
    ring_sum.assign(nbins, 0);
@@ -20,6 +20,7 @@ AdaptiveSpotFinderCPU::AdaptiveSpotFinderCPU(const AzimuthalIntegrationMapping &
    ring_sigma.assign(nbins, 0.0f);
    ring_thr.assign(nbins, 0.0f);
    ring_bkg.assign(nbins, NAN);
+    ring_bits.assign(OutputSize(), 0);
 }

 // Accumulate per-ring mean/variance from the raw (photon) image. clip_k <= 0 -> use every valid
@@ -61,9 +62,7 @@ void AdaptiveSpotFinderCPU::AccumulateRings(const ImagePreprocessorBuffer &image

 void AdaptiveSpotFinderCPU::Detect(const ImagePreprocessorBuffer &image,
                                   const SpotFindingSettings &settings) {
-    const auto &pixel_to_bin = mapping.GetPixelToBin();
    const size_t nbins = ring_sum.size();
-    const size_t npix = static_cast<size_t>(width) * height;

    // --- Stage A: robust per-ring background (one plain pass + two sigma-clip passes) ---
    AccumulateRings(image, 0.0f);
@@ -104,9 +103,28 @@ void AdaptiveSpotFinderCPU::Detect(const ImagePreprocessorBuffer &image,
                          ? g_thr
                          : adaptive_threshold::RingThreshold(ring_mean[b], ring_sigma[b], p, z);

-    // --- Stage C: flag strong pixels into the bit buffer (value >= ring threshold) ---
+    // --- Stage C: the ring threshold, intersected with the classic local-box SNR test ---
+    FlagRings(image);
+
+    if (settings.signal_to_noise_threshold <= 0.0f) {
+        // No local test asked for: the ring threshold alone decides, as the fixed photon floor
+        // alone would in the classic finder.
+        output_buffer = ring_bits;
+        return;
+    }
+
+    // The ring threshold IS the photon floor here, so the local pass must not apply another one.
+    SpotFindingSettings local = settings;
+    local.photon_count_threshold = 0;
+    ImageSpotFinderCPU::Detect(image, local);
    for (size_t i = 0; i < OutputSize(); ++i)
-        output_buffer[i] = 0;
+        output_buffer[i] &= ring_bits[i];
+}
+
+void AdaptiveSpotFinderCPU::FlagRings(const ImagePreprocessorBuffer &image) {
+    const auto &pixel_to_bin = mapping.GetPixelToBin();
+    const size_t nbins = ring_thr.size();
+    const size_t npix = static_cast<size_t>(width) * height;

    std::bitset<32> out = 0;
    for (size_t pxl = 0; pxl < npix; ++pxl) {
@@ -122,10 +140,10 @@ void AdaptiveSpotFinderCPU::Detect(const ImagePreprocessorBuffer &image,
        if (strong)
            out.set(bit);
        if (bit == 31) {
-            output_buffer[pxl / 32] = out.to_ulong();
+            ring_bits[pxl / 32] = out.to_ulong();
            out.reset();
        }
    }
    if (npix % 32 != 0)
-        output_buffer[OutputSize() - 1] = out.to_ulong();
+        ring_bits[OutputSize() - 1] = out.to_ulong();
 }
@@ -5,7 +5,7 @@

 #include <vector>

-#include "ImageSpotFinder.h"
+#include "ImageSpotFinderCPU.h"
 #include "SpotFindingSettings.h"
 #include "../../common/AzimuthalIntegrationMapping.h"

@@ -22,9 +22,29 @@
 // whose (peak-excluded) background mean is mu, the threshold is the smallest count whose Poisson
 // upper tail is <= p = E / N_pixels, max'd with a Gaussian arm mu + z*sigma to absorb read/flat-field
 // excess. Because it is set from the image's own noise, the SAME E lands at ~12 photons on the first
-// set and ~5 on the weaker one with no user input. Detection then is simply value > ring_threshold,
-// fed to the same connected-component builder as the classic finder.
-class AdaptiveSpotFinderCPU : public ImageSpotFinder {
+// set and ~5 on the weaker one with no user input.
+//
+// The ring threshold replaces the floor and ONLY the floor: the classic finder's local-box SNR test
+// still has to pass, which is why this engine runs it (ImageSpotFinderCPU) and intersects the two
+// masks.
+//
+// A whole-ring threshold is an ABSOLUTE contour with no feedback from the pixel's own surroundings,
+// so the area a spot puts above it grows as sigma^2 * ln(peak/threshold) and never saturates: on a
+// strongly diffracting rotation set the detected footprint grows by 8 pixels per e-fold of peak, so
+// the brightest reflections came out as 100-500 pixel blobs. The local box has no such contour. The
+// spot inflates the box's own variance, and the peak divides out of the acceptance test, so the box
+// cuts every spot at roughly a fixed FRACTION of its own height - a peak-relative contour. Measured
+// on the same set, the footprint then grows by -0.2 pixels per e-fold, i.e. not at all, and lands on
+// the classic finder's own number to two decimals.
+//
+// The two arms bind in different regimes, which is the point of intersecting rather than choosing.
+// On serial stills the ring background is a fraction of a count and the ring threshold lands BELOW
+// the fixed floor the classic finder would use, so the ring arm decides and the local box passes
+// everything - which is the whole reason this engine exists. On a bright rotation set the ring
+// background is tens of counts, the ring threshold lands several times ABOVE that floor, and the
+// local box decides. The engine is therefore never worse than the classic finder on footprint, and
+// never worse than a fixed floor on a weak background.
+class AdaptiveSpotFinderCPU : public ImageSpotFinderCPU {
    const AzimuthalIntegrationMapping &mapping;

    // per-ring scratch, sized to the mapping's bin count
@@ -40,8 +60,13 @@ class AdaptiveSpotFinderCPU : public ImageSpotFinder {
    // ring_mean of the last Detect(), NaN where the ring holds too few pixels to be its own background.
    // Kept separately because ring_mean carries the previous frame's value for an empty ring.
    std::vector<float> ring_bkg;
+    // Pixels at or above their ring's threshold, packed like output_buffer. Intersected with the
+    // local-box mask that ImageSpotFinderCPU::Detect leaves in output_buffer.
+    std::vector<uint32_t> ring_bits;

    void AccumulateRings(const ImagePreprocessorBuffer &image, float clip_k);
+    // Fill ring_bits from the thresholds of the current frame.
+    void FlagRings(const ImagePreprocessorBuffer &image);

 public:
    explicit AdaptiveSpotFinderCPU(const AzimuthalIntegrationMapping &mapping);
@@ -276,14 +276,20 @@ __global__ void flag_strong(const int32_t *__restrict__ image,
    }
 }

+// out &= mask, over the packed bit buffer. One word per thread; the buffer is one bit per pixel, so
+// this is 1/32 of an image pass.
+__global__ void and_bits(uint32_t *__restrict__ out, const uint32_t *__restrict__ mask, size_t nwords) {
+    for (size_t i = blockIdx.x * blockDim.x + threadIdx.x; i < nwords; i += blockDim.x * gridDim.x)
+        out[i] &= mask[i];
+}
+
 } // namespace

 AdaptiveSpotFinderGPU::AdaptiveSpotFinderGPU(const AzimuthalIntegrationMapping &in_mapping,
                                             std::shared_ptr<CudaStream> in_stream)
-    : ImageSpotFinder(static_cast<int32_t>(in_mapping.GetWidth()),
-                      static_cast<int32_t>(in_mapping.GetHeight()), false),
+    : ImageSpotFinderGPU(static_cast<int32_t>(in_mapping.GetWidth()),
+                         static_cast<int32_t>(in_mapping.GetHeight()), std::move(in_stream)),
      mapping(in_mapping),
-      stream(in_stream),
      nbins(in_mapping.GetBinNumber()),
      npix(in_mapping.GetPixelToBin().size()),
      gpu_sum(nbins),
@@ -294,7 +300,7 @@ AdaptiveSpotFinderGPU::AdaptiveSpotFinderGPU(const AzimuthalIntegrationMapping &
      gpu_sum_corr(nbins),
      gpu_sum2_corr(nbins),
      gpu_thr(nbins),
-      gpu_strong(OutputSize()),
+      gpu_ring(OutputSize()),
      host_sum(nbins),
      host_sum2(nbins),
      host_count(nbins),
@@ -308,8 +314,6 @@ AdaptiveSpotFinderGPU::AdaptiveSpotFinderGPU(const AzimuthalIntegrationMapping &
      prof_sum_reg(prof_sum),
      prof_sum2_reg(prof_sum2),
      prof_count_reg(prof_count),
-      extractor(static_cast<int32_t>(in_mapping.GetWidth()),
-                static_cast<int32_t>(in_mapping.GetHeight()), std::move(in_stream)),
      last_profile(in_mapping) {

    // The current device, not device 0: callers round-robin engines across GPUs, so device 0's shared
@@ -464,27 +468,30 @@ void AdaptiveSpotFinderGPU::Detect(const ImagePreprocessorBuffer &image,

    if (host_thr.empty()) {
        // Nothing valid to threshold against: leave no strong pixels for the extractor to build on.
-        cuda_err(cudaMemsetAsync(gpu_strong, 0, OutputByteSize(), *stream));
+        cuda_err(cudaMemsetAsync(gpu_out_1, 0, OutputByteSize(), *stream));
        cuda_err(cudaStreamSynchronize(*stream));
        return;
    }

-    // --- Stage C: flag strong pixels into the bit buffer (value >= ring threshold) ---
+    // --- Stage C: the ring threshold, intersected with the classic local-box SNR test ---
    cuda_err(cudaMemcpyAsync(gpu_thr, host_thr.data(), sizeof(float) * nbins, cudaMemcpyHostToDevice, *stream));
-    cuda_err(cudaMemsetAsync(gpu_strong, 0, OutputByteSize(), *stream));
+    cuda_err(cudaMemsetAsync(gpu_ring, 0, OutputByteSize(), *stream));
    flag_strong<<<flag_blocks, flag_threads, 0, *stream>>>(
-        image.getGPUBuffer(), gpu_pixel_to_bin->get(), gpu_thr, gpu_strong, npix, nbins);
+        image.getGPUBuffer(), gpu_pixel_to_bin->get(), gpu_thr, gpu_ring, npix, nbins);
+
+    if (settings.signal_to_noise_threshold <= 0.0f) {
+        // No local test asked for: the ring threshold alone decides, as the fixed photon floor
+        // alone would in the classic finder.
+        cuda_err(cudaMemcpyAsync(gpu_out_1, gpu_ring, OutputByteSize(), cudaMemcpyDeviceToDevice, *stream));
+        cuda_err(cudaStreamSynchronize(*stream));
+        return;
+    }
+
+    // The ring threshold IS the photon floor here, so the local pass must not apply another one.
+    SpotFindingSettings local = settings;
+    local.photon_count_threshold = 0;
+    ImageSpotFinderGPU::Detect(image, local);
+    and_bits<<<flag_blocks, flag_threads, 0, *stream>>>(gpu_out_1, gpu_ring, OutputSize());
    // The bit buffer stays on the device and ExtractComponents reads it there, on this same stream,
-    // so the ordering already guarantees flag_strong has finished. Waiting here only idled the host.
-}
-
-void AdaptiveSpotFinderGPU::SetResolutionMaskBits(const std::vector<uint32_t> &packed_mask) {
-    ImageSpotFinder::SetResolutionMaskBits(packed_mask);
-    extractor.SetResolutionMask(res_mask_bits);
-}
-
-const std::vector<DiffractionSpot> &AdaptiveSpotFinderGPU::ExtractComponents(const ImagePreprocessorBuffer &image,
-                                                                             const SpotFindingSettings &settings) {
-    extractor.Extract(gpu_strong, image.getGPUBuffer(), settings, components);
-    return components;
+    // so the ordering already guarantees and_bits has finished. Waiting here only idled the host.
 }
@@ -12,20 +12,23 @@
 //   - the azimuthal-integration profile (mean intensity per ring, in flat-field-corrected space), and
 //   - the per-ring background (mean, sigma, peak-excluded via two sigma-clip passes) that sets the
 //     self-calibrating spot-detection threshold (in raw photon counts).
-// It then flags strong pixels (value >= ring threshold) into a packed bit buffer and hands that
-// buffer - still on the device - to SpotExtractorGPU, which builds the spots there.
+// The ring threshold is the photon-count FLOOR and nothing more: this engine derives from the classic
+// ImageSpotFinderGPU and intersects the ring mask with that engine's local-box SNR mask, exactly as
+// AdaptiveSpotFinderCPU does on the host (the reason the local test is not optional is written out
+// there). The result is left in the inherited bit buffer, still on the device, and the inherited
+// SpotExtractorGPU builds the spots from it.
 //
 // Numerically it reproduces AdaptiveSpotFinderCPU: the same three-pass robust background, the same
 // per-ring threshold formula (shared via AdaptiveThreshold.h, computed on the host once per frame),
-// and the same raw-count detection test. The only differences from the CPU are those inherent to a
-// GPU reduction (float per-ring accumulation in atomic order vs the CPU's serial double sums), which
-// shift a handful of borderline pixels at most. The corrected sums for the azint profile are
+// and the same intersection with the local test. The only differences from the CPU are those inherent
+// to a GPU reduction (float per-ring accumulation in atomic order vs the CPU's serial double sums),
+// which shift a handful of borderline pixels at most. The corrected sums for the azint profile are
 // accumulated in the SAME plain first pass, so one reduction feeds both products.

 #include <memory>
 #include <vector>

-#include "ImageSpotFinder.h"
+#include "ImageSpotFinderGPU.h"
 #include "SpotExtractorGPU.h"
 #include "SpotFindingSettings.h"
 #include "../../common/AzimuthalIntegrationProfile.h"
@@ -33,9 +36,8 @@
 #include "../indexing/CUDAMemHelpers.h"
 #include "../indexing/CudaSharedTables.h"

-class AdaptiveSpotFinderGPU : public ImageSpotFinder {
+class AdaptiveSpotFinderGPU : public ImageSpotFinderGPU {
    const AzimuthalIntegrationMapping &mapping;
-    std::shared_ptr<CudaStream> stream;

    const int nbins;
    const size_t npix;
@@ -68,9 +70,10 @@ class AdaptiveSpotFinderGPU : public ImageSpotFinder {
    CudaDevicePtr<float>    gpu_sum_corr;
    CudaDevicePtr<float>    gpu_sum2_corr;

-    // Per-ring detection threshold (host-computed, uploaded) and the strong-pixel bit buffer.
+    // Per-ring detection threshold (host-computed, uploaded) and the ring-threshold mask that is
+    // intersected into the inherited bit buffer.
    CudaDevicePtr<float>    gpu_thr;
-    CudaDevicePtr<uint32_t> gpu_strong;
+    CudaDevicePtr<uint32_t> gpu_ring;

    // Host mirrors of the small per-ring transfers.
    std::vector<unsigned long long> host_sum;    // clipped raw sum   } input to the host threshold computation
@@ -94,8 +97,6 @@ class AdaptiveSpotFinderGPU : public ImageSpotFinder {
    CudaRegisteredVector<float>    prof_sum2_reg;
    CudaRegisteredVector<uint32_t> prof_count_reg;

-    SpotExtractorGPU extractor;   // builds the spots from gpu_strong without it leaving the device
-
    AzimuthalIntegrationProfile last_profile;   // filled every Run(), retrievable via GetProfile()

    // One reduction pass over the image into the raw accumulators. clip_k <= 0 -> plain pass (all
@@ -114,10 +115,6 @@ public:
    AdaptiveSpotFinderGPU &operator=(const AdaptiveSpotFinderGPU &) = delete;

    void Detect(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) override;
-    void SetResolutionMaskBits(const std::vector<uint32_t> &packed_mask) override;
-    [[nodiscard]] uint32_t StrongPixelCount() const override { return extractor.StrongPixelCount(); }
-    const std::vector<DiffractionSpot> &ExtractComponents(const ImagePreprocessorBuffer &image,
-                                                          const SpotFindingSettings &settings) override;

    // The azimuthal profile computed as a byproduct of the last Detect() - lets this engine replace the
    // separate azint pass in the analysis pipeline.
@@ -61,10 +61,27 @@ inline double NormalQuantile(double p) {
 // significance floor while the background is countable (it carries the sqrt(mu) shot-noise
 // implicitly, so a bright low-resolution ring gets a high threshold). It DEGENERATES at mu -> 0
 // (a single photon on a zero background is "significant"), which is why it is max'd with a
-// read-noise-floored Gaussian arm by the caller. Short-circuits to Gaussian for large mu.
+// read-noise-floored Gaussian arm by the caller.
+//
+// Past the summation limit the quantile is taken from the Cornish-Fisher expansion
+// (Cornish and Fisher (1938) Rev. Int. Stat. Inst. 5, 307-320), whose skewness
+// term (z^2-1)/6 is what a plain mu + z*sqrt(mu) leaves out. At the 4-6 sigma this operating point
+// works at, that term is 3-6 counts, so the Gaussian form alone stood BELOW the true Poisson
+// quantile - and it did so with a step at the switch, since below it the exact quantile was used.
+// Cornish-Fisher is within one count of the exact value at every mu, so the two arms now join
+// smoothly.
+//
+// How much this is worth depends on which arm of RingThreshold wins, and on measured data it is
+// often neither: where the ring background is over-dispersed (a clipped ring sigma that still
+// carries the ring's own azimuthal structure, 1.2-4.9x sqrt(mu) on a strongly diffracting rotation
+// set) the Gaussian arm is the larger of the two on every ring above mu = 50 and this correction
+// changes no threshold at all. It is the right value to return regardless: a caller that ever sees
+// the Poisson arm win up there would otherwise get a bar that jumps at mu = 50.
 inline float PoissonThreshold(double mu, double p, double z) {
-    if (mu > 50.0)
-        return static_cast<float>(mu + z * std::sqrt(mu));
+    // The summation below needs k up to about mu + z*sqrt(mu), and exp(-mu) has to stay normal.
+    constexpr double SUM_LIMIT = 200.0;
+    if (mu > SUM_LIMIT)
+        return static_cast<float>(mu + z * std::sqrt(mu) + (z * z - 1.0) / 6.0);
    if (mu < 1e-6) mu = 1e-6;
    const double target = 1.0 - p;
    double pmf = std::exp(-mu);
@@ -11,12 +11,18 @@
 #include "../indexing/CUDAMemHelpers.h"

 class ImageSpotFinderGPU : public ImageSpotFinder {
+protected:
+    // Protected rather than private because AdaptiveSpotFinderGPU derives from this engine: it is
+    // this same local-box detection with the fixed photon floor replaced by a per-resolution-ring
+    // one, so it reuses the stream, the bit buffers and the extractor rather than owning a second
+    // set of them.
    std::shared_ptr<CudaStream> stream;

    CudaDevicePtr<uint32_t> gpu_out_0;
    CudaDevicePtr<uint32_t> gpu_out_1;   // holds the strong-pixel bits after Detect()
    SpotExtractorGPU extractor;

+private:
    const int numberOfCudaThreads = 128;    // #threads per block that should work well for Nvidia L4
    const int numberOfWaves = 32;           // #waves that should work well for Nvidia L4
    const int windowSizeLimit = 32;         // limit on the window size (2nby+1, 2nbx+1) to prevent shared memory problems
@@ -32,10 +32,19 @@ constexpr int64_t SPOT_SHAPE_FREE_PIXELS = 50;
 //
 // It is INERT on every dataset it has been measured on: gate on and gate off give a byte-identical
 // merge on three rotation crystals including the strongly diffracting one the raised bound was written
-// for. What it guards against is a population none of them carry in quantity - the raised bound admits
-// components up to 200 px, and on ice-arc or cosmic-track-heavy data those are what arrives. So it is
-// kept as a bound on the SHAPE of what the larger size bound now lets through, not because any
-// measurement here needed it.
+// for. The raised bound is what did nearly all of the work there - of the components an absolute ring
+// contour pushes past 50 pixels, this shape test rejects a few percent and the raise re-admits the
+// rest.
+//
+// Both only matter while a spot's footprint can grow with its brightness, and since the adaptive
+// finder intersects its ring threshold with the local-box SNR test the footprint is peak-relative and
+// no component on that set reaches 50 pixels at all.
+//
+// The shape test is kept anyway because it is the one of the two that does not go stale: a size bound
+// is a bet on how large spots are, and the detectors and the detection rule both move underneath it -
+// under every peak-relative detector tried, the 200-pixel bound became unreachable while this test
+// still fired. The bound is kept alongside it because the cap the local box imposes scales with the
+// spot WIDTH, so a set with wider spots than the ones measured here will reach past 50 pixels again.
 // A shape test rather than a size test is what XDS and DISTL bound spots with; neither uses this
 // statistic, but the choice of shape over size is theirs.
 // Following Kabsch (2010) Acta Cryst. D66, 125-132 and Zhang, Sauter et al. (2006) J. Appl. Cryst. 39, 112-119
@@ -1,6 +1,7 @@
 // SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
 // SPDX-License-Identifier: GPL-3.0-only

+#include <cmath>
 #include <algorithm>

 #include <catch2/catch_all.hpp>
@@ -190,3 +191,64 @@ TEST_CASE("AdaptiveSpotFinderCPU_RingBackgroundExcludesPeaks", "[AdaptiveSpotFin
    for (size_t b = 0; b < clean.size(); b++)
        CHECK(std::isfinite(clean[b]) == (pixels_per_bin[b] >= 40));
 }
+
+// The footprint of a bright reflection must not grow with its brightness. A whole-ring threshold is
+// an absolute contour, so the area a Gaussian puts above it grows as sigma^2 ln(peak/threshold) - the
+// same spot detected a hundred times brighter comes back tens of pixels larger, and an upper bound on
+// spot size becomes an upper bound on spot INTENSITY. Intersecting with the local-box test removes
+// that: the spot inflates the box's own variance, the peak divides out of the acceptance test, and the
+// contour lands at a fixed fraction of the spot's own height whatever that height is.
+TEST_CASE("AdaptiveSpotFinderCPU_FootprintDoesNotGrowWithBrightness", "[AdaptiveSpotFinder]") {
+    DiffractionExperiment x(DetJF4M());
+    x.DetectorDistance_mm(80).BeamX_pxl(1030).BeamY_pxl(1080);
+    x.QSpacingForAzimInt_recipA(0.05).QRangeForAzimInt_recipA(0.05, 5.0);
+    x.GeometryTransformation(false);
+
+    PixelMask pixel_mask(x);
+    AzimuthalIntegrationMapping mapping(x, pixel_mask);
+    const auto &pixel_to_bin = mapping.GetPixelToBin();
+
+    const size_t w = x.GetXPixelsNum();
+    const size_t h = x.GetYPixelsNum();
+
+    size_t spot_row = 0, spot_col = 0;
+    for (size_t row = 400; row < h - 400 && spot_row == 0; row++)
+        for (size_t col = 400; col < w - 400; col++)
+            if (pixel_to_bin[row * w + col] != UINT16_MAX) {
+                spot_row = row;
+                spot_col = col;
+                break;
+            }
+    REQUIRE(spot_row > 0);
+
+    std::vector<bool> res_mask(x.GetPixelsNum(), false);
+    auto settings = AdaptiveSettings();
+    settings.min_pix_per_spot = 2;
+    settings.max_pix_per_spot = 100000;   // no bound, so the footprint itself is what is measured
+
+    // One Gaussian of width 1.5 px on a flat background of 10, at three amplitudes a hundred apart.
+    auto footprint = [&](double amplitude) {
+        ImagePreprocessorBuffer buffer(x.GetPixelsNum());
+        for (size_t i = 0; i < w * h; i++)
+            buffer[i] = 10;
+        constexpr double sigma = 1.5;
+        for (int dr = -12; dr <= 12; dr++)
+            for (int dc = -12; dc <= 12; dc++) {
+                const double r2 = dr * dr + dc * dc;
+                buffer[(spot_row + dr) * w + spot_col + dc] =
+                    10 + static_cast<int32_t>(amplitude * std::exp(-r2 / (2 * sigma * sigma)));
+            }
+        AdaptiveSpotFinderCPU finder(mapping);
+        finder.SetResolutionMask(res_mask);
+        const auto spots = finder.Run(buffer, settings);
+        REQUIRE(spots.size() == 1);
+        return spots[0].PixelCount();
+    };
+
+    const int64_t small = footprint(300.0);
+    const int64_t large = footprint(30000.0);
+    CHECK(small > 0);
+    // A hundredfold in peak is 4.6 e-folds. An absolute contour would add sigma^2 ln(100) ~ 10 pixels
+    // per e-fold of AREA here, tens of pixels in all; a peak-relative one adds nothing.
+    CHECK(large - small <= 4);
+}
@@ -51,7 +51,7 @@ TEST_CASE("AdaptiveThreshold_PoissonThreshold", "[SpotFinding]") {
    const float z = static_cast<float>(NormalQuantile(1.0 - p));

    // The defining property: the returned count is the SMALLEST whose upper tail is within p.
-    for (const double mu: {1e-6, 0.1, 1.0, 3.0, 10.0, 40.0}) {
+    for (const double mu: {1e-6, 0.1, 1.0, 3.0, 10.0, 40.0, 60.0, 120.0}) {
        const int thr = static_cast<int>(PoissonThreshold(mu, p, z));
        CHECK(PoissonUpperTail(mu, thr) <= p);
        CHECK(PoissonUpperTail(mu, thr - 1) > p);
@@ -65,8 +65,19 @@ TEST_CASE("AdaptiveThreshold_PoissonThreshold", "[SpotFinding]") {
        prev = thr;
    }

-    // Above mu = 50 it short-circuits to the Gaussian form mu + z sqrt(mu).
-    CHECK(PoissonThreshold(100.0, p, z) == Catch::Approx(100.0 + z * 10.0).epsilon(1e-5));
+    // Above the summation limit the Cornish-Fisher form takes over, and it has to stay a POISSON
+    // quantile: the skewness term (z^2-1)/6 is what a plain mu + z*sqrt(mu) leaves out, and at this
+    // z the Gaussian form alone lets through several times the tail asked for.
+    for (const double mu: {250.0, 400.0}) {
+        const float thr = PoissonThreshold(mu, p, z);
+        CHECK(thr > mu + z * std::sqrt(mu));
+        CHECK(PoissonUpperTail(mu, static_cast<int>(thr)) <= 2 * p);
+        CHECK(PoissonUpperTail(mu, static_cast<int>(mu + z * std::sqrt(mu))) > 2 * p);
+    }
+
+    // And it joins the exact quantile smoothly at the switch - no step for a ring whose background
+    // drifts across it from frame to frame.
+    CHECK(std::fabs(PoissonThreshold(200.5, p, z) - PoissonThreshold(199.5, p, z)) < 2.0f);

    // A tighter operating point (smaller p) can only raise the threshold.
    CHECK(PoissonThreshold(5.0, 1e-8, static_cast<float>(NormalQuantile(1.0 - 1e-8)))