diff --git a/acquisition_device/AcquisitionDevice.cpp b/acquisition_device/AcquisitionDevice.cpp index 9c2e7160..ed410ab9 100644 --- a/acquisition_device/AcquisitionDevice.cpp +++ b/acquisition_device/AcquisitionDevice.cpp @@ -324,6 +324,6 @@ void AcquisitionDevice::SetSpotFinderParameters(const SpotFindingSettings &setti fpga_parameters.count_threshold = settings.photon_count_threshold; fpga_parameters.max_d = settings.low_resolution_limit; fpga_parameters.min_d = settings.high_resolution_limit; - fpga_parameters.min_pix_per_spot = settings.min_pix_per_spot; + fpga_parameters.min_pix_per_spot = settings.min_pix_per_spot.value_or(2); HW_SetSpotFinderParameters(fpga_parameters); } diff --git a/broker/OpenAPIConvert.cpp b/broker/OpenAPIConvert.cpp index 626fd716..4783dcb2 100644 --- a/broker/OpenAPIConvert.cpp +++ b/broker/OpenAPIConvert.cpp @@ -32,7 +32,7 @@ org::openapitools::server::model::Spot_finding_settings Convert(const SpotFindin org::openapitools::server::model::Spot_finding_settings ret; ret.setSignalToNoiseThreshold(input.signal_to_noise_threshold); ret.setPhotonCountThreshold(input.photon_count_threshold); - ret.setMinPixPerSpot(input.min_pix_per_spot); + ret.setMinPixPerSpot(input.min_pix_per_spot.value_or(2)); ret.setMaxPixPerSpot(input.max_pix_per_spot); ret.setHighResolutionLimit(input.high_resolution_limit); ret.setLowResolutionLimit(input.low_resolution_limit); diff --git a/common/DiffractionExperiment.cpp b/common/DiffractionExperiment.cpp index 81593874..38f50880 100644 --- a/common/DiffractionExperiment.cpp +++ b/common/DiffractionExperiment.cpp @@ -608,8 +608,8 @@ void DiffractionExperiment::CheckDataProcessingSettings(const SpotFindingSetting check_finite("Signal to noise threshold", settings.signal_to_noise_threshold); check_min("Signal to noise threshold", settings.signal_to_noise_threshold, 1); check_min("Photon count threshold", settings.photon_count_threshold, 0); - check_min("Minimum pixels per spot", settings.min_pix_per_spot, 1); - check_min("Maximum pixels per spot", settings.max_pix_per_spot, settings.min_pix_per_spot + 1); + check_min("Minimum pixels per spot", settings.min_pix_per_spot.value_or(2), 1); + check_min("Maximum pixels per spot", settings.max_pix_per_spot, settings.min_pix_per_spot.value_or(2) + 1); check_finite("Spot finding high resolution limit", settings.high_resolution_limit); check_finite("Spot finding low resolution limit", settings.low_resolution_limit); diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index d9aa7df5..c284ac71 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -219,6 +219,14 @@ Strong pixels are grouped into connected components (adjacent strong pixels) usi Spot-level filters include minimum/maximum pixel count and resolution limits. +### 3.5 Adaptive per-image minimum spot size + +The minimum-pixels-per-spot filter (§3.4) trades sensitivity against noise: a small value keeps faint one- or two-pixel spots — real signal on strong data, but detector noise on high-background frames — while a larger value keeps only well-formed spots. The best value is dataset-dependent, so for serial-stills indexing it can be chosen **per image** rather than fixed. The frame is indexed three times, at min-pix 3, 2 and 1, and the setting that maximises + +$$ \frac{n_\mathrm{indexed}^2}{n_\mathrm{total}} \quad\text{(indexed-spot count weighted by indexed fraction)} $$ + +is kept; the frame is then integrated once at that min-pix. The fraction factor discounts the extra spots a smaller min-pix admits *unless the lattice actually explains them*, so strong frames keep their real weak spots (extending resolution) while noise-flooded frames stay strict. Because min-pix filters the connected components *after* detection, the three attempts only repeat the cheap CCL and spot-level filter, not the pixel reduction; the azimuthal profile is identical across them. This is a **stills-only, indexing-path** option — rotation indexing builds one global lattice from all frames and keeps a fixed min-pix. In `rugnux` it is the default; giving an explicit `--min-pix-per-spot` pins a fixed value instead. + --- ## 4. Indexing overview diff --git a/docs/RUGNUX.md b/docs/RUGNUX.md index 953c5034..68075a24 100644 --- a/docs/RUGNUX.md +++ b/docs/RUGNUX.md @@ -150,15 +150,18 @@ then merge against a reference structure: rugnux serial_master.h5 \ -o serial_run -N 32 \ -X ffbidx -C 79,79,38,90,90,90 -S 96 \ - --spot-sigma 4 \ + --adaptive-spots \ -z reference.mtz \ --scaling-high-resolution 1.8 ``` `ffbidx` requires a known cell (`-C`) and is the indexer of choice for sparse serial stills. For -weak serial data, tightening spot finding with `--spot-sigma 4` typically raises the indexing rate -substantially. If a dataset *does* carry a goniometer axis but you want per-frame stills processing -anyway, add `--force-still`. +serial stills, prefer the self-calibrating spot finder (`--adaptive-spots`) and leave +`--min-pix-per-spot` **unset** so it is chosen per image — across the still-target battery this +combination raises the indexing rate and typically extends resolution over a fixed threshold and +fixed min-pix, at equal or better CC½. (You can still pin a fixed threshold with `--spot-sigma` / +`--spot-threshold` and a fixed min-pix with `--min-pix-per-spot`.) If a dataset *does* carry a +goniometer axis but you want per-frame stills processing anyway, add `--force-still`. ## Command-line options @@ -186,9 +189,11 @@ Spot finding: | --- | --- | | `--spot-sigma ` | Noise sigma level for spot finding (default: 3.0) | | `--spot-threshold ` | Photon-count threshold for spot finding (default: 10) | +| `--adaptive-spots` | Self-calibrating detection: replace the fixed `--spot-threshold` with a per-resolution-ring threshold derived from each image's own noise, so one setting adapts across datasets (no per-dataset `--spot-threshold`/`--spot-sigma` tuning) | +| `--spot-false-pixels ` | Adaptive-detection operating point: expected noise pixels tolerated per frame (default: 100; implies `--adaptive-spots`) | | `--spot-high-resolution ` | High-resolution limit for spot finding, Å (default: 1.5) | | `--spot-low-resolution ` | Low-resolution limit for spot finding, Å (default: 50; lower it, e.g. 24, to exclude the direct-beam halo on weak serial data) | -| `--min-pix-per-spot ` | Minimum connected strong pixels per spot (default: 2; serial data can index better with 1 and a higher `--spot-threshold`) | +| `--min-pix-per-spot ` | Minimum connected strong pixels per spot. **If omitted, min-pix is chosen per image** (stills indexing): the frame is indexed at min-pix 3/2/1 and the one maximising indexed-spot count × indexed fraction is kept. Give an explicit value to force a fixed min-pix instead. | | `--max-spots ` | Maximum spot count (default: 250) | | `--detect-ice-rings[=on\|off]` | Flag ice-ring spots (de-prioritised in indexing) and exclude ice-ring reflections from scaling/merging; overrides the dataset/master-file setting (default: use the dataset value) | diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 9b503036..08956bec 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -510,13 +510,14 @@ IndexAndRefine::DetermineRefineAnalyze(DataMessage &msg, const SpotFindingSettin return outcome; } -void IndexAndRefine::ProcessImage(DataMessage &msg, +bool IndexAndRefine::ProcessImage(DataMessage &msg, const SpotFindingSettings &spot_finding_settings, BraggPrediction &prediction, const BraggIntegrateFn &integrate) { auto outcome = DetermineRefineAnalyze(msg, spot_finding_settings); if (outcome && spot_finding_settings.quick_integration) QuickPredictAndIntegrate(msg, spot_finding_settings, prediction, integrate, *outcome); + return outcome.has_value(); } bool IndexAndRefine::IndexFrameOnly(DataMessage &msg, const SpotFindingSettings &spot_finding_settings) { diff --git a/image_analysis/IndexAndRefine.h b/image_analysis/IndexAndRefine.h index e870c6b2..dda9a0d3 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -102,7 +102,9 @@ public: prediction_mosaicity_override_ = std::move(mosaicity_per_frame); } - void ProcessImage(DataMessage &msg, const SpotFindingSettings &settings, + // Returns whether the frame indexed (a lattice was found and refined). Integration, when it runs, + // is a further step gated on quick_integration. + bool ProcessImage(DataMessage &msg, const SpotFindingSettings &settings, BraggPrediction &prediction, const BraggIntegrateFn &integrate); // Index a single frame (no integration) with the current forced rotation lattice; used to score // first-pass sampling schemes on the real per-image path. Returns whether the frame indexed. diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index bf021e89..7d188a38 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -3,6 +3,8 @@ #include "MXAnalysisWithoutFPGA.h" +#include + #include "spot_finding/StrongPixelSet.h" #include "../compression/JFJochDecompress.h" @@ -117,30 +119,68 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, || mask_low_res != spot_finding_settings.low_resolution_limit) UpdateMaskResolution(spot_finding_settings); - const auto spot_finding_start_time = std::chrono::steady_clock::now(); ImageSpotFinder &finder = spot_finding_settings.adaptive_threshold ? static_cast(*adaptiveSpotFinder) : *spotFinder; - const std::vector spots = finder.Run(*preprocessor_buffer, spot_finding_settings, mask_resolution); - SpotAnalyze(experiment, spot_finding_settings, spots, output); - const auto spot_finding_end_time = std::chrono::steady_clock::now(); - output.spot_finding_time_s = std::chrono::duration(spot_finding_end_time - spot_finding_start_time).count(); + const auto integrate_fn = [this](const std::vector &predicted, size_t npredicted, + int64_t image_number) { + return bragg_engine->Run(*preprocessor_buffer, predicted, npredicted, image_number); + }; + + // A missing min-pix (std::nullopt) means "choose it per image". This applies only to the stills + // indexing path (each frame is indexed independently); rotation indexing builds one lattice from + // all frames, so it keeps the fixed min-pix and the single-pass finder. + const bool adaptive_min_pix = !spot_finding_settings.min_pix_per_spot.has_value() + && spot_finding_settings.indexing + && !experiment.IsRotationIndexing(); + if (adaptive_min_pix) { + // Choose the per-image min-pix adaptively instead of a fixed one. min-pix filters + // connected components AFTER detection, so re-running the finder only re-does the cheap CCL + + // spot filter, not the reduction; the azimuthal profile is identical across attempts. Index + // at 3/2/1 (index-only, no integration/accumulation) and keep whichever maximises + // n_indexed^2 / n_total (indexed count weighted by indexed fraction), then integrate once at + // that min-pix. spot_finding_time_s covers the whole escalation. + const auto start_time = std::chrono::steady_clock::now(); + SpotFindingSettings s = spot_finding_settings; + int best_mp = 0; + double best_score = -1.0; + for (int mp : {3, 2, 1}) { + s.min_pix_per_spot = mp; + const std::vector spots = finder.Run(*preprocessor_buffer, s, mask_resolution); + SpotAnalyze(experiment, s, spots, output); + if (indexer.IndexFrameOnly(output, s)) { + const double n_idx = static_cast(output.spot_count_indexed.value_or(0)); + const double n_tot = static_cast(std::max(1, output.spot_count.value_or(1))); + const double score = n_idx * n_idx / n_tot; + if (score > best_score) { best_score = score; best_mp = mp; } + } + } + if (best_mp != 0) { + // Re-run spot finding + index at the winning min-pix and integrate there. + s.min_pix_per_spot = best_mp; + const std::vector spots = finder.Run(*preprocessor_buffer, s, mask_resolution); + SpotAnalyze(experiment, s, spots, output); + indexer.ProcessImage(output, s, *prediction, integrate_fn); + } + output.spot_finding_time_s = std::chrono::duration(std::chrono::steady_clock::now() - start_time).count(); + } else { + const auto spot_finding_start_time = std::chrono::steady_clock::now(); + const std::vector spots = finder.Run(*preprocessor_buffer, spot_finding_settings, mask_resolution); + SpotAnalyze(experiment, spot_finding_settings, spots, output); + output.spot_finding_time_s = std::chrono::duration(std::chrono::steady_clock::now() - spot_finding_start_time).count(); + if (spot_finding_settings.indexing) + indexer.ProcessImage(output, spot_finding_settings, *prediction, integrate_fn); + } #ifdef JFJOCH_USE_CUDA if (fused) { - // Lift the azimuthal profile the fused engine computed in the same pass; its azint cost is - // folded into spot_finding_time_s above. + // Lift the azimuthal profile the fused engine computed in the same pass (identical across + // any min-pix retries); its azint cost is folded into spot_finding_time_s above. profile.Clear(integration); profile += fused_adaptive->GetProfile(); output.azint_time_s = 0.0f; } #endif - - if (spot_finding_settings.indexing) - indexer.ProcessImage(output, spot_finding_settings, *prediction, - [this](const std::vector &predicted, size_t npredicted, int64_t image_number) { - return bragg_engine->Run(*preprocessor_buffer, predicted, npredicted, image_number); - }); } output.max_viable_pixel_value = ret.max_value; diff --git a/image_analysis/spot_finding/SpotFindingSettings.h b/image_analysis/spot_finding/SpotFindingSettings.h index b69b7072..96a4c9af 100644 --- a/image_analysis/spot_finding/SpotFindingSettings.h +++ b/image_analysis/spot_finding/SpotFindingSettings.h @@ -10,7 +10,11 @@ struct SpotFindingSettings { bool enable = true; float signal_to_noise_threshold = 4.0; // STRONG_PIXEL in XDS int64_t photon_count_threshold = 10; // Threshold in photon counts - int64_t min_pix_per_spot = 2; // Minimum pixels per spot + // Minimum connected pixels per spot. std::nullopt = choose it per image: on the stills indexing + // path the frame is indexed at min-pix 3/2/1 and the one maximising indexed count x indexed fraction + // 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 float high_resolution_limit = 2.0; float low_resolution_limit = 50.0; diff --git a/image_analysis/spot_finding/StrongPixelSet.cpp b/image_analysis/spot_finding/StrongPixelSet.cpp index 439d3050..3a6ac524 100644 --- a/image_analysis/spot_finding/StrongPixelSet.cpp +++ b/image_analysis/spot_finding/StrongPixelSet.cpp @@ -89,7 +89,7 @@ void StrongPixelSet::FindSpotsImage(const SpotFindingSettings &settings, std::ve if (!pixels.empty() && (strong_pixel_count < UINT16_MAX)) { for (const auto &spot: sparseccl()) { if ((spot.PixelCount() <= settings.max_pix_per_spot) - && (spot.PixelCount() >= settings.min_pix_per_spot)) { + && (spot.PixelCount() >= settings.min_pix_per_spot.value_or(2))) { spots.push_back(spot); } } @@ -102,7 +102,7 @@ void StrongPixelSet::FindSpots(const DiffractionExperiment &experiment, const Sp if (!pixels.empty() && (strong_pixel_count < UINT16_MAX)) { for (const auto &spot: sparseccl()) { if ((spot.PixelCount() <= settings.max_pix_per_spot) - && (spot.PixelCount() >= settings.min_pix_per_spot)) { + && (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/rugnux/RugnuxCommandLine.cpp b/rugnux/RugnuxCommandLine.cpp index e2fdab15..751821fe 100644 --- a/rugnux/RugnuxCommandLine.cpp +++ b/rugnux/RugnuxCommandLine.cpp @@ -87,6 +87,10 @@ std::string RugnuxCommandLine(const ProcessConfig &config, add("--spot-threshold", std::to_string(sf.photon_count_threshold)); if (sf.adaptive_threshold) add("--spot-false-pixels", num(sf.false_pixels_per_frame)); + // min-pix is chosen per image unless an explicit value is given, so emit --min-pix-per-spot only + // when a fixed min-pix was selected; its absence selects the adaptive per-image path. + if (sf.min_pix_per_spot.has_value()) + add("--min-pix-per-spot", std::to_string(*sf.min_pix_per_spot)); add("--spot-high-resolution", num(sf.high_resolution_limit)); add("--max-spots", std::to_string(experiment.GetMaxSpotCount())); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 17eff6e5..15b579d7 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -6,6 +6,7 @@ #include #include #include +#include #include #include #include @@ -64,7 +65,7 @@ void print_usage() { std::cout << " Spot finding" << std::endl; std::cout << " --spot-sigma Noise sigma level for spot finding (default: 3.0)" << std::endl; std::cout << " --spot-threshold Photon count threshold for spot finding (default: 10)" << std::endl; - std::cout << " --min-pix-per-spot Minimum connected strong pixels per spot (default: 2; serial data can index better with 1 + a higher --spot-threshold)" << std::endl; + std::cout << " --min-pix-per-spot Minimum connected strong pixels per spot. If omitted, min-pix is chosen PER IMAGE (stills indexing): the frame is indexed at min-pix 3/2/1 and the one maximising indexed count x indexed fraction is kept. Give an explicit value to force a fixed min-pix instead." << std::endl; std::cout << " --adaptive-spots Self-calibrating detection: replace the fixed --spot-threshold with a per-resolution-ring threshold set from each image's own noise, so one setting adapts across datasets (no per-dataset --spot-threshold/--spot-sigma tuning)" << std::endl; std::cout << " --spot-false-pixels Adaptive detection operating point: expected noise pixels tolerated per frame (default: 100; implies --adaptive-spots)" << std::endl; std::cout << " --spot-high-resolution High resolution limit for spot finding (default: 1.5)" << std::endl; @@ -513,7 +514,7 @@ int main(int argc, char **argv) { std::optional max_spot_count_override; float sigma_spot_finding = 3.0; int64_t photon_count_threshold_spot_finding = 10; - int64_t min_pix_per_spot = 2; + std::optional min_pix_per_spot; // unset -> adaptive per image; a value -> fixed min-pix bool adaptive_spots = false; float false_pixels_per_frame = 100.0f; bool refine_bfactor = false; @@ -741,8 +742,9 @@ int main(int argc, char **argv) { photon_count_threshold_spot_finding); break; case OPT_MIN_PIX_PER_SPOT: + // Giving an explicit min-pix opts out of the per-image adaptive selection. min_pix_per_spot = parse_number_arg(optarg, "--min-pix-per-spot", logger, 1); - logger.Info("Minimum pixels per spot set to {:d}", min_pix_per_spot); + logger.Info("Minimum pixels per spot fixed at {:d} (adaptive per-image min-pix off)", *min_pix_per_spot); break; case OPT_ADAPTIVE_SPOTS: adaptive_spots = true; diff --git a/viewer/widgets/JFJochViewerSettingsDock.cpp b/viewer/widgets/JFJochViewerSettingsDock.cpp index 82379ae0..3203b55b 100644 --- a/viewer/widgets/JFJochViewerSettingsDock.cpp +++ b/viewer/widgets/JFJochViewerSettingsDock.cpp @@ -231,12 +231,18 @@ QWidget *JFJochViewerSettingsDock::BuildMXPage() { "settings are not used while this is on."); auto *highResSpot = new SliderPlusBox(0.5, 5.0, 0.1, 1, page); highResSpot->setValue(spot_.high_resolution_limit); - auto *minPix = new NumberLineEdit(1.0f, 50.0f, static_cast(spot_.min_pix_per_spot), 0, "px", page); + auto *minPix = new NumberLineEdit(1.0f, 50.0f, static_cast(spot_.min_pix_per_spot.value_or(2)), 0, "px", page); + auto *adaptiveMinPix = new QCheckBox("Adaptive min-pix (per image)", page); + adaptiveMinPix->setChecked(!spot_.min_pix_per_spot.has_value()); + adaptiveMinPix->setToolTip("Choose the minimum pixels/spot per image (stills indexing): index at " + "min-pix 3/2/1 and keep whichever maximises indexed count x indexed " + "fraction. The fixed min-pixels/spot value is not used while this is on."); auto *maxSpots = new NumberLineEdit(10.0f, 100000.0f, static_cast(max_spots_), 0, "", page); spot->addRow("", adaptive); spot->addRow("Signal/noise", snr); spot->addRow("Photon count", count); spot->addRow("High resolution [Å]", highResSpot); + spot->addRow("", adaptiveMinPix); spot->addRow("Min pixels/spot", minPix); spot->addRow("Max spots/image", maxSpots); spotSection->setContentLayout(spot); @@ -264,6 +270,18 @@ QWidget *JFJochViewerSettingsDock::BuildMXPage() { syncAdaptiveEnabled(on); EmitSpotFinding(); }); + // Adaptive min-pix chooses the value per image (min_pix_per_spot = std::nullopt), so the fixed + // min-pixels/spot field is unused while it is on - grey it out to make that clear. + auto syncMinPixEnabled = [minPix](bool adaptive_on) { minPix->setEnabled(!adaptive_on); }; + syncMinPixEnabled(!spot_.min_pix_per_spot.has_value()); + connect(adaptiveMinPix, &QCheckBox::toggled, this, [this, minPix, syncMinPixEnabled](bool on) { + if (on) + spot_.min_pix_per_spot = std::nullopt; + else + spot_.min_pix_per_spot = std::llround(minPix->value()); + syncMinPixEnabled(on); + EmitSpotFinding(); + }); // --- Indexing --- auto *idxSection = new CollapsibleSection("Indexing", page);