Files
Jungfraujoch/tests/StrongPixelSetTest.cpp
T
leonarski_fandClaude Opus 5 0af109b838
Build Packages / build:rugnux:aarch64 (cross) (push) Successful in 8m39s
Build Packages / build:windows:nocuda (push) Successful in 17m3s
Build Packages / build:windows:cuda (push) Successful in 19m14s
Build Packages / build:rugnux-tgz (x86_64) (push) Successful in 19m38s
Build Packages / build:viewer-tgz:cpu (push) Successful in 20m32s
Build Packages / build:rpm (rocky9_nocuda) (push) Successful in 23m30s
Build Packages / build:viewer-tgz:cuda (push) Successful in 23m43s
Build Packages / build:rpm (ubuntu2204_nocuda) (push) Successful in 27m8s
Build Packages / build:rpm (rocky8_nocuda) (push) Successful in 27m12s
Build Packages / build:rugnux:windows (push) Successful in 10m52s
Build Packages / build:rpm (ubuntu2404_nocuda) (push) Successful in 19m27s
Build Packages / build:rpm (rocky9_sls9) (push) Successful in 21m13s
Build Packages / build:rpm (rocky8_sls9) (push) Successful in 26m49s
Build Packages / build:rpm (rocky9) (push) Successful in 23m36s
Build Packages / build:rpm (ubuntu2404) (push) Successful in 24m9s
Build Packages / Generate python client (push) Successful in 32s
Build Packages / build:rpm (rocky8) (push) Successful in 29m17s
Build Packages / Create release (push) Skipped
Build Packages / XDS test (durin plugin) (push) Successful in 11m18s
Build Packages / Build documentation (push) Successful in 1m10s
Build Packages / DIALS test (push) Successful in 25m11s
Build Packages / build:rpm (ubuntu2204) (push) Successful in 28m6s
Build Packages / XDS test (neggia plugin) (push) Successful in 8m56s
Build Packages / XDS test (JFJoch plugin) (push) Successful in 9m51s
Build Packages / Unit tests (push) Successful in 1h26m0s
Spot finding: bound how large a spot may be, not how bright
A component of more than 50 connected pixels was discarded. Under the
self-calibrating threshold that is an INTENSITY CEILING, not a size bound: the
threshold is an absolute per-ring contour, so a component's area above it grows
as sigma^2*ln(A/T), without bound in the peak amplitude. Measured on the strong
rotation set, footprints run 3 px at 30-100 counts to 50 px above 10000 - a
slope of 4.4 px per ln(peak) against the fixed local-box test's 0.8, which
saturates at 13 px and never reaches the bound at all. So the brighter a
reflection, the more certainly it was thrown away: 59% of the box finder's
d>3 A spots were missing from the adaptive finder's list, including every one
of its ten strongest, at an intensity ratio of 1.085 for those that did match.
Indexed spots per image collapsed from 220 to 7.

Raise the bound to 200 - CrystFEL peakfinder8's --max-pix-count, the only
directly comparable number in the field; XDS has no such parameter and guards
on shape instead - and ask a component above 50 pixels to be COMPACT: it must
fill 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, and those are what
an upper bound was ever protecting against. Below 50 nothing is asked of the
shape, so every component accepted before still is. Integer arithmetic on both
sides, and on the GPU the bounding side fits in what was padding, so the device
struct does not grow.

The shape test is what makes the raise safe. With a flat 200 alone, two battery
crystals moved: one lost a little I/sigma, and the other's de-novo lattice was
NOT MONOTONE in the bound - correct at 50, 100 and 200, wrong at 150 and at 250
and above - so 200 was partly luck. With the shape test both are unchanged to
three decimals.

De novo at defaults on the strong set: indexing rate 0.574 -> 0.804,
completeness 97.3 -> 99.6%, <I/sigma> 3.42 -> 6.07, R_meas 0.299 -> 0.249,
CC1/2 0.947 -> 0.961, ISa 3.29 -> 4.13. The full 37-crystal battery, scored per
shell, is still owed.

Also documents what StrongPixelSet::AddStrongPixel has required since the
component search became linear - pixels in raster order - and puts the existing
test's insertion order into it. Both callers scan a bitmap in ascending flat
index and always satisfied it; the test did not, and was the only thing that
did not.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01FBumeJVx4oeXxiBRpkrE5H
2026-08-28 11:49:49 +02:00

155 lines
6.0 KiB
C++

// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include <catch2/catch_all.hpp>
#include "../image_analysis/spot_finding/StrongPixelSet.h"
TEST_CASE("DiffractionSpot_AddOperator","[StrongPixelSet]") {
DiffractionSpot spot1(4,1,10), spot2(3,4,5);
spot1 += spot2;
// Spot2 is not changed
REQUIRE(spot2.Count() == 5);
REQUIRE(spot1.Count() == 15);
REQUIRE(spot1.RawCoord().x == Catch::Approx((4*10+3*5)/15.0));
REQUIRE(spot1.RawCoord().y == Catch::Approx((1*10+4*5)/15.0));
REQUIRE(spot1.MaxCount() == 10);
REQUIRE(spot1.PixelCount() == 2);
}
TEST_CASE("StrongPixelSet_BuildSpots","[StrongPixelSet]") {
DiffractionExperiment experiment(DetJF(1,1,0,0,false));
experiment.Raw();
SpotFindingSettings settings;
settings.low_resolution_limit = 200.0;
settings.high_resolution_limit = 0.5;
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<DiffractionSpot> spots;
StrongPixelSet strong_pixel_set;
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);
REQUIRE(spots[0].Count() == 9.0f);
REQUIRE(spots[0].PixelCount() == 9);
REQUIRE(spots[0].RawCoord().x == Catch::Approx(7.0));
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<DiffractionSpot> 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));
experiment.BeamX_pxl(1000).BeamY_pxl(1000).DetectorDistance_mm(100).PhotonEnergy_keV(WVL_1A_IN_KEV);
StrongPixelSet strong_pixel_set(experiment);
DeviceOutput output{};
for (auto &i: output.spot_finding_result.strong_pixel)
i = 0;
output.pixels[0] = 7;
output.pixels[1] = 7;
output.pixels[1024] = 7;
output.pixels[1025] = 7;
output.pixels[7] = 345;
// Pixels 0, 1, 7, 1024 and 1025 will be enabled
output.spot_finding_result.strong_pixel[0] = (1U<<0) | (1U<<1) | (1U<<7);
output.spot_finding_result.strong_pixel[RAW_MODULE_COLS/8] = (1<<0) | (1<<1);
output.spot_finding_result.strong_pixel_count = 5;
strong_pixel_set.ReadFPGAOutput(output);
REQUIRE(strong_pixel_set.GetStrongPixelCount() == 5);
SpotFindingSettings settings{
.min_pix_per_spot = 2,
.max_pix_per_spot = 55,
.high_resolution_limit = 1.0,
.low_resolution_limit = 50.0
};
std::vector<DiffractionSpot> spots;
strong_pixel_set.FindSpots(experiment, settings, spots, 2);
CHECK(spots.size() == 1);
REQUIRE(spots[0].RawCoord().x == Approx(0.5));
REQUIRE(spots[0].RawCoord().y == Approx(CONVERTED_MODULE_LINES * 3 + 36 * 2 - 1 - 0.5));
REQUIRE(spots[0].PixelCount() == 4);
REQUIRE(spots[0].MaxCount() == 7);
}
TEST_CASE("StrongPixelSet_ReadFPGAOutput_2","[StrongPixelSet]") {
DiffractionExperiment experiment(DetectorGeometry(8, 2, 8, 36, true));
experiment.BeamX_pxl(1000).BeamY_pxl(1000).DetectorDistance_mm(100).PhotonEnergy_keV(WVL_1A_IN_KEV);
StrongPixelSet strong_pixel_set(experiment);
DeviceOutput output{};
for (auto &i: output.spot_finding_result.strong_pixel)
i = 0;
output.pixels[0] = 7;
output.pixels[1] = 7;
output.pixels[1024] = 7;
output.pixels[1025] = 7;
output.pixels[7] = 345;
// Pixels 0, 1, 7, 1024 and 1025 will be enabled
output.spot_finding_result.strong_pixel[0] = (1U<<0) | (1U<<1) | (1U<<7);
output.spot_finding_result.strong_pixel[RAW_MODULE_COLS/8] = (1<<0) | (1<<1);
output.spot_finding_result.strong_pixel_count = 5;
strong_pixel_set.ReadFPGAOutput(output);
REQUIRE(strong_pixel_set.GetStrongPixelCount() == 5);
SpotFindingSettings settings{
.min_pix_per_spot = 1,
.max_pix_per_spot = 3,
.high_resolution_limit = 1.0,
.low_resolution_limit = 50.0
};
std::vector<DiffractionSpot> spots;
strong_pixel_set.FindSpots(experiment, settings, spots, 2);
CHECK(spots.size() == 1);
REQUIRE(spots[0].RawCoord().x == Approx(7));
REQUIRE(spots[0].RawCoord().y == Approx(CONVERTED_MODULE_LINES * 3 + 36 * 2 - 1));
REQUIRE(spots[0].PixelCount() == 1);
REQUIRE(spots[0].MaxCount() == 345);
} */