diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index b2940fb9c..37030a2aa 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -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. diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp index a60e78327..9f7cd16ad 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp +++ b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp @@ -9,8 +9,8 @@ #include "AdaptiveThreshold.h" AdaptiveSpotFinderCPU::AdaptiveSpotFinderCPU(const AzimuthalIntegrationMapping &in_mapping) - : ImageSpotFinder(static_cast(in_mapping.GetWidth()), - static_cast(in_mapping.GetHeight())), + : ImageSpotFinderCPU(static_cast(in_mapping.GetWidth()), + static_cast(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(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(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(); } diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h index c78d732b0..3604d1b9c 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h +++ b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h @@ -5,7 +5,7 @@ #include -#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 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 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); diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu index 65ed5123b..2cafc5c06 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu +++ b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu @@ -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 in_stream) - : ImageSpotFinder(static_cast(in_mapping.GetWidth()), - static_cast(in_mapping.GetHeight()), false), + : ImageSpotFinderGPU(static_cast(in_mapping.GetWidth()), + static_cast(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(in_mapping.GetWidth()), - static_cast(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<<>>( - 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<<>>(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 &packed_mask) { - ImageSpotFinder::SetResolutionMaskBits(packed_mask); - extractor.SetResolutionMask(res_mask_bits); -} - -const std::vector &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. } diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h index c15d48c0e..0a3bcb545 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h +++ b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h @@ -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 #include -#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 stream; const int nbins; const size_t npix; @@ -68,9 +70,10 @@ class AdaptiveSpotFinderGPU : public ImageSpotFinder { CudaDevicePtr gpu_sum_corr; CudaDevicePtr 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 gpu_thr; - CudaDevicePtr gpu_strong; + CudaDevicePtr gpu_ring; // Host mirrors of the small per-ring transfers. std::vector host_sum; // clipped raw sum } input to the host threshold computation @@ -94,8 +97,6 @@ class AdaptiveSpotFinderGPU : public ImageSpotFinder { CudaRegisteredVector prof_sum2_reg; CudaRegisteredVector 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 &packed_mask) override; - [[nodiscard]] uint32_t StrongPixelCount() const override { return extractor.StrongPixelCount(); } - const std::vector &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. diff --git a/image_analysis/spot_finding/AdaptiveThreshold.h b/image_analysis/spot_finding/AdaptiveThreshold.h index 047ee4698..745c64d27 100644 --- a/image_analysis/spot_finding/AdaptiveThreshold.h +++ b/image_analysis/spot_finding/AdaptiveThreshold.h @@ -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(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(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); diff --git a/image_analysis/spot_finding/ImageSpotFinderGPU.h b/image_analysis/spot_finding/ImageSpotFinderGPU.h index d63f4acd1..1d9382771 100644 --- a/image_analysis/spot_finding/ImageSpotFinderGPU.h +++ b/image_analysis/spot_finding/ImageSpotFinderGPU.h @@ -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 stream; CudaDevicePtr gpu_out_0; CudaDevicePtr 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 diff --git a/image_analysis/spot_finding/SpotFindingSettings.h b/image_analysis/spot_finding/SpotFindingSettings.h index 8b79aab6c..897408920 100644 --- a/image_analysis/spot_finding/SpotFindingSettings.h +++ b/image_analysis/spot_finding/SpotFindingSettings.h @@ -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 diff --git a/tests/AdaptiveSpotFinderCPUTest.cpp b/tests/AdaptiveSpotFinderCPUTest.cpp index 819ccbf1b..0db0270a6 100644 --- a/tests/AdaptiveSpotFinderCPUTest.cpp +++ b/tests/AdaptiveSpotFinderCPUTest.cpp @@ -1,6 +1,7 @@ // SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute // SPDX-License-Identifier: GPL-3.0-only +#include #include #include @@ -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 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(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); +} diff --git a/tests/AdaptiveThresholdTest.cpp b/tests/AdaptiveThresholdTest.cpp index ecd13da64..dc5c88257 100644 --- a/tests/AdaptiveThresholdTest.cpp +++ b/tests/AdaptiveThresholdTest.cpp @@ -51,7 +51,7 @@ TEST_CASE("AdaptiveThreshold_PoissonThreshold", "[SpotFinding]") { const float z = static_cast(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(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(thr)) <= 2 * p); + CHECK(PoissonUpperTail(mu, static_cast(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(NormalQuantile(1.0 - 1e-8)))