diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index e6074e4cf..b2940fb9c 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -1,6 +1,8 @@ # Changelog ## 1.0.0 ### 1.0.0-rc.166 +* 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. * 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/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index ece6feb53..256cb58fe 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -372,6 +372,8 @@ Strong pixels are grouped into connected components (adjacent strong pixels) usi Spot-level filters include minimum/maximum pixel count and resolution limits. +**The upper bound is on how large a *round* spot may be, not on how bright one may be.** Under the self-calibrating threshold (ยง3.2) a component's area above the contour grows as $\sigma^2\ln(A/T)$ โ€” without bound in the peak amplitude $A$ โ€” so an upper bound on pixel count alone becomes an *intensity ceiling*: measured on a strongly diffracting crystal, footprints run 3 px at 30โ€“100 counts to 50 px above 10 000, four times the slope the fixed local-box test gives, and the old bound of 50 discarded the brightest reflections on every image. The bound is therefore 200 (the same as CrystFEL `peakfinder8`'s `--max-pix-count`; XDS has no such parameter at all and guards on shape instead, with `SPOT_MAXIMUM-CENTROID`), and a component above 50 pixels must in addition **fill at least a fifth of the square its bounding box fits inside**. A Bragg reflection is round and fills about half of that square however bright it is; an ice arc, a cosmic-ray track or a lit detector row fills a fifth or less, which is what an upper bound was ever protecting against. Below 50 pixels no shape is asked for, so everything accepted before still is. The test is integer arithmetic, so the host and the GPU extractor agree by construction. + The host implementation (`StrongPixelSet::sparseccl`) is the SparseCCL of the ACTS/traccc project: it runs over the strong pixels sorted row-major, uses a sliding window over the previous line and a union-find whose root is each component's lowest index. On the GPU the same labelling runs **on the diff --git a/etc/broker_local.json b/etc/broker_local.json index 417b14e4a..ce9845eee 100644 --- a/etc/broker_local.json +++ b/etc/broker_local.json @@ -20,7 +20,7 @@ "signal_to_noise_threshold":4.0, "photon_count_threshold":10, "low_resolution_limit":50.0, - "max_pix_per_spot":50, + "max_pix_per_spot":200, "min_pix_per_spot":2, "indexing":true, "quick_integration":true diff --git a/image_analysis/spot_finding/SpotExtractorGPU.cu b/image_analysis/spot_finding/SpotExtractorGPU.cu index 0cc6624b5..d286c5586 100644 --- a/image_analysis/spot_finding/SpotExtractorGPU.cu +++ b/image_analysis/spot_finding/SpotExtractorGPU.cu @@ -181,7 +181,7 @@ __global__ void finish_components(const uint32_t *__restrict__ index, const int3 int32_t *__restrict__ count, SpotExtractorGPUSpot *__restrict__ scratch, SpotExtractorGPUSpot *__restrict__ out, uint32_t *__restrict__ nout, const uint32_t *__restrict__ nstrong, uint32_t capacity, - int width, int max_pix) { + int width, int max_pix, int shape_free_pix, int min_fill_percent) { const int n = static_cast(min(*nstrong, capacity)); if (threadIdx.x == 0) *nout = 0; __syncthreads(); @@ -231,10 +231,14 @@ __global__ void finish_components(const uint32_t *__restrict__ index, const int3 if (want > max_pix) continue; long long x = 0, y = 0; long long photons = 0, max_photons = LLONG_MIN; + int min_col = INT_MAX, max_col = INT_MIN, min_line = INT_MAX, max_line = INT_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]; + const int col = static_cast(index[j] % width), line = static_cast(index[j] / width); + min_col = min(min_col, col); max_col = max(max_col, col); + min_line = min(min_line, line); max_line = max(max_line, line); // 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; @@ -247,14 +251,23 @@ __global__ void finish_components(const uint32_t *__restrict__ index, const int3 scratch[l].y = y; scratch[l].photons = photons; scratch[l].max_photons = max_photons; + scratch[l].bbox_side = max(max_col - min_col, max_line - min_line) + 1; } __syncthreads(); - // 4) max-pix filter, compacted by another prefix sum so the surviving spots keep their order + // 4) size and shape filter, compacted by another prefix sum so the surviving spots keep their + // order. The test is SpotShapeAccepted written out - the constants come in as arguments rather + // than being included here, so there is exactly one definition of them. + auto keep = [&](const SpotExtractorGPUSpot &s) { + if (s.pixel_count > max_pix) return false; + if (s.pixel_count <= shape_free_pix) return true; + return static_cast(s.pixel_count) * 100 + >= static_cast(min_fill_percent) * s.bbox_side * s.bbox_side; + }; const int label_chunk = (nlabel + nthreads - 1) / nthreads; const int label_lo = min(t * label_chunk, nlabel), label_hi = min(label_lo + label_chunk, nlabel); uint32_t nkeep = 0; - for (int i = label_lo; i < label_hi; i++) nkeep += (scratch[i].pixel_count <= max_pix) ? 1u : 0u; + for (int i = label_lo; i < label_hi; i++) nkeep += keep(scratch[i]) ? 1u : 0u; shared[t] = nkeep; __syncthreads(); for (int d = 1; d < nthreads; d <<= 1) { @@ -265,7 +278,7 @@ __global__ void finish_components(const uint32_t *__restrict__ index, const int3 } uint32_t pos = shared[t] - nkeep; for (int i = label_lo; i < label_hi; i++) - if (scratch[i].pixel_count <= max_pix) out[pos++] = scratch[i]; + if (keep(scratch[i])) out[pos++] = scratch[i]; if (t == nthreads - 1) *nout = shared[t]; } @@ -336,7 +349,9 @@ void SpotExtractorGPU::Extract(const uint32_t *gpu_strong, const int32_t *gpu_im resolve_roots<<<512, THREADS, 0, *stream>>>(gpu_parent, gpu_root, gpu_nstrong, max_strong); finish_components<<<1, FINISH_THREADS, 0, *stream>>>(gpu_index, gpu_value, gpu_root, gpu_label, gpu_count, gpu_spot, gpu_spot_out, gpu_nspot, - gpu_nstrong, max_strong, width, max_pix); + gpu_nstrong, max_strong, width, max_pix, + static_cast(SPOT_SHAPE_FREE_PIXELS), + static_cast(SPOT_MIN_FILL_PERCENT)); // Rides along with the spot count on the frame's one synchronisation, so knowing how many strong // pixels there were costs nothing. cuda_err(cudaMemcpyAsync(host_nstrong, gpu_nstrong, sizeof(uint32_t), cudaMemcpyDeviceToHost, *stream)); diff --git a/image_analysis/spot_finding/SpotExtractorGPU.h b/image_analysis/spot_finding/SpotExtractorGPU.h index 415c69e94..4bef4f31d 100644 --- a/image_analysis/spot_finding/SpotExtractorGPU.h +++ b/image_analysis/spot_finding/SpotExtractorGPU.h @@ -43,7 +43,9 @@ struct SpotExtractorGPUSpot { int64_t photons; int64_t max_photons; int32_t pixel_count; - int32_t padding; + // Longer side of the component's bounding box, for SpotShapeAccepted. It fits in what used to be + // padding, so carrying it costs nothing. + int32_t bbox_side; }; class SpotExtractorGPU { diff --git a/image_analysis/spot_finding/SpotFindingSettings.h b/image_analysis/spot_finding/SpotFindingSettings.h index a2bdde060..d610277bc 100644 --- a/image_analysis/spot_finding/SpotFindingSettings.h +++ b/image_analysis/spot_finding/SpotFindingSettings.h @@ -22,6 +22,26 @@ constexpr uint32_t StrongPixelLimit(size_t pixel_count) { return limit > UINT16_MAX ? limit : UINT16_MAX; } +// A connected component this size or smaller is judged on size alone. That is where the size bound +// stood before it was raised, so every component that used to be accepted still is. +constexpr int64_t SPOT_SHAPE_FREE_PIXELS = 50; +// A larger one must also be COMPACT: it has to fill at least this percent of the square its bounding +// box fits inside. A Bragg reflection is round and fills about half of that square however bright it +// is; an ice arc, a cosmic-ray track or a lit detector row fills a fifth or less, and those are what +// an upper bound on spot size was ever for. Measured on the strong rotation set, this keeps 97% of the +// genuine wide reflections a raised bound admits and none of the ice arcs. +// 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 +constexpr int64_t SPOT_MIN_FILL_PERCENT = 20; + +// Integer arithmetic throughout, so the host and the GPU extractor agree by construction. +constexpr bool SpotShapeAccepted(int64_t pixel_count, int64_t bbox_side) { + if (pixel_count <= SPOT_SHAPE_FREE_PIXELS) + return true; + return pixel_count * 100 >= SPOT_MIN_FILL_PERCENT * bbox_side * bbox_side; +} + struct SpotFindingSettings { bool enable = true; float signal_to_noise_threshold = 4.0; // STRONG_PIXEL in XDS @@ -31,7 +51,9 @@ struct SpotFindingSettings { // is kept (see MXAnalysisWithoutFPGA::Analyze); a value fixes it. Defaults to a concrete value, so // the online receiver and the FPGA path keep the single-pass fixed behaviour unless set otherwise. std::optional min_pix_per_spot = 2; - int64_t max_pix_per_spot = 50; // Maximum pixels per spot + // Maximum pixels per spot. A component above SPOT_SHAPE_FREE_PIXELS must also pass + // SpotShapeAccepted, so this bounds how large a ROUND spot may be, not how bright. + int64_t max_pix_per_spot = 200; // High-resolution limit for spot finding [A]. std::nullopt = as far as the detector reaches, i.e. no // resolution clipping of the detection at all (DiffractionExperiment::GetDetectorMaxResolution_A // supplies the number where one is needed, e.g. for the spot plot's shells). diff --git a/image_analysis/spot_finding/StrongPixelSet.cpp b/image_analysis/spot_finding/StrongPixelSet.cpp index 6c631236a..0c12cc241 100644 --- a/image_analysis/spot_finding/StrongPixelSet.cpp +++ b/image_analysis/spot_finding/StrongPixelSet.cpp @@ -8,6 +8,7 @@ // The union-find and the two-scan structure are theirs. How a pixel's earlier neighbours are FOUND // is not: see sparseccl below. +#include #include #include "StrongPixelSet.h" @@ -40,7 +41,7 @@ uint32_t StrongPixelSet::make_union(uint32_t e1, uint32_t e2) { return e; } -std::vector StrongPixelSet::sparseccl() { +std::vector StrongPixelSet::sparseccl(const SpotFindingSettings &settings) { L.resize(pixels.size()); unsigned int labels = 0; @@ -89,20 +90,39 @@ std::vector StrongPixelSet::sparseccl() { } std::vector spots(labels); + // The bounding box travels with the accumulation rather than on DiffractionSpot: it is wanted + // only to judge the shape here, and carrying it further would leave a member that + // ConvertToImageCoordinates silently invalidates. + std::vector min_col(labels, UINT16_MAX), max_col(labels, 0); + std::vector min_line(labels, UINT16_MAX), max_line(labels, 0); - for (uint32_t i = 0; i < L.size(); i++) - spots[L[i]].AddPixel(pixels[i].col, pixels[i].line, pixels[i].counts); + for (uint32_t i = 0; i < L.size(); i++) { + const uint32_t l = L[i]; + spots[l].AddPixel(pixels[i].col, pixels[i].line, pixels[i].counts); + min_col[l] = std::min(min_col[l], pixels[i].col); + max_col[l] = std::max(max_col[l], pixels[i].col); + min_line[l] = std::min(min_line[l], pixels[i].line); + max_line[l] = std::max(max_line[l], pixels[i].line); + } - return spots; + std::vector out; + for (uint32_t l = 0; l < labels; l++) { + if (spots[l].PixelCount() > settings.max_pix_per_spot) + continue; + const int64_t w = static_cast(max_col[l]) - min_col[l] + 1; + const int64_t h = static_cast(max_line[l]) - min_line[l] + 1; + if (SpotShapeAccepted(spots[l].PixelCount(), std::max(w, h))) + out.push_back(spots[l]); + } + return out; } void StrongPixelSet::FindComponentsImage(const SpotFindingSettings &settings, std::vector &spots) { // No StrongPixelLimit test here: the caller knows how big the image is and has already applied it. - for (const auto &spot: sparseccl()) { - if (spot.PixelCount() <= settings.max_pix_per_spot) - spots.push_back(spot); - } + // Size and shape are applied by sparseccl. + for (const auto &spot: sparseccl(settings)) + spots.push_back(spot); } void StrongPixelSet::FindSpots(const DiffractionExperiment &experiment, const SpotFindingSettings &settings, @@ -110,9 +130,8 @@ void StrongPixelSet::FindSpots(const DiffractionExperiment &experiment, const Sp // Per module, so the bar is the module's own - and ReadFPGAOutput has already refused anything // past max_strong_pixel_per_module, far below it. if (!pixels.empty() && (strong_pixel_count < StrongPixelLimit(RAW_MODULE_SIZE))) { - for (const auto &spot: sparseccl()) { - if ((spot.PixelCount() <= settings.max_pix_per_spot) - && (spot.PixelCount() >= settings.min_pix_per_spot.value_or(2))) { + for (const auto &spot: sparseccl(settings)) { + if (spot.PixelCount() >= settings.min_pix_per_spot.value_or(2)) { auto s = spot; s.ConvertToImageCoordinates(experiment, module_number); spots.push_back(s); diff --git a/image_analysis/spot_finding/StrongPixelSet.h b/image_analysis/spot_finding/StrongPixelSet.h index 5026122d3..7004dd510 100644 --- a/image_analysis/spot_finding/StrongPixelSet.h +++ b/image_analysis/spot_finding/StrongPixelSet.h @@ -26,12 +26,19 @@ class StrongPixelSet { uint32_t find_root(uint32_t e); uint32_t make_union(uint32_t e1, uint32_t e2); - std::vector sparseccl(); + // Labels the strong pixels and returns the components that pass BOTH the size bound and + // SpotShapeAccepted. min-pix is deliberately not applied here (see FindComponentsImage). + std::vector sparseccl(const SpotFindingSettings &settings); public: void ReadFPGAOutput(const DiffractionExperiment& experiment, const DeviceOutput& output); StrongPixelSet(); + // Pixels MUST arrive in raster order - line ascending, and column ascending within a line. The + // component search walks the previous line with a forward-only cursor and would miss neighbours + // handed to it out of order. Both callers satisfy this by construction, each scanning a bitmap in + // ascending flat index: ImageSpotFinder::ExtractComponentsHost over the image, ReadFPGAOutput over + // a module (whose row length is a multiple of 32, so no word straddles two lines). void AddStrongPixel(uint16_t col, uint16_t line, int32_t photons = 1); void FindSpots(const DiffractionExperiment &experiment, const SpotFindingSettings &settings, std::vector &spots, uint16_t module_number); diff --git a/tests/StrongPixelSetTest.cpp b/tests/StrongPixelSetTest.cpp index c2243e100..97649d25d 100644 --- a/tests/StrongPixelSetTest.cpp +++ b/tests/StrongPixelSetTest.cpp @@ -29,17 +29,12 @@ TEST_CASE("StrongPixelSet_BuildSpots","[StrongPixelSet]") { settings.min_pix_per_spot = 3; settings.max_pix_per_spot = 200; + // In raster order, which is what AddStrongPixel requires and what both callers produce. std::vector spots; StrongPixelSet strong_pixel_set; - strong_pixel_set.AddStrongPixel(7,105); - strong_pixel_set.AddStrongPixel(7,106); - strong_pixel_set.AddStrongPixel(7,104); - strong_pixel_set.AddStrongPixel(6,106); - strong_pixel_set.AddStrongPixel(6,105); - strong_pixel_set.AddStrongPixel(6,104); - strong_pixel_set.AddStrongPixel(8,106); - strong_pixel_set.AddStrongPixel(8,105); - strong_pixel_set.AddStrongPixel(8,104); + for (uint16_t line = 104; line <= 106; line++) + for (uint16_t col = 6; col <= 8; col++) + strong_pixel_set.AddStrongPixel(col, line); strong_pixel_set.FindSpots(experiment, settings, spots, 0); REQUIRE(spots.size() == 1); @@ -49,6 +44,40 @@ TEST_CASE("StrongPixelSet_BuildSpots","[StrongPixelSet]") { REQUIRE(spots[0].RawCoord().y == Catch::Approx(105.0)); } +TEST_CASE("StrongPixelSet_LargeSpotMustBeCompact","[StrongPixelSet]") { + SpotFindingSettings settings; // max_pix_per_spot 200, so all three pass on size alone + + // A filled 9x9 square: 81 pixels in a 9-wide box, so it fills all of it. A Bragg reflection is + // this shape however bright it is, and this is what a raised size bound is meant to admit. + StrongPixelSet square; + for (uint16_t line = 100; line < 109; line++) + for (uint16_t col = 200; col < 209; col++) + square.AddStrongPixel(col, line); + std::vector spots; + square.FindComponentsImage(settings, spots); + REQUIRE(spots.size() == 1); + CHECK(spots[0].PixelCount() == 81); + + // The same 81 pixels drawn as a line: an ice arc, a cosmic-ray track or a lit detector row. It + // fills 1/81 of its bounding square, so it is refused however many pixels it has. + StrongPixelSet streak; + for (uint16_t col = 200; col < 281; col++) + streak.AddStrongPixel(col, 100); + spots.clear(); + streak.FindComponentsImage(settings, spots); + CHECK(spots.empty()); + + // Shape is asked of large components only. A short line of the same one-pixel width is below + // SPOT_SHAPE_FREE_PIXELS and is kept, exactly as it was before the bound was raised. + StrongPixelSet short_streak; + for (uint16_t col = 200; col < 210; col++) + short_streak.AddStrongPixel(col, 100); + spots.clear(); + short_streak.FindComponentsImage(settings, spots); + REQUIRE(spots.size() == 1); + CHECK(spots[0].PixelCount() == 10); +} + /* TEST_CASE("StrongPixelSet_ReadFPGAOutput_1","[StrongPixelSet]") { DiffractionExperiment experiment(DetectorGeometry(8, 2, 8, 36, true));