From 16bf3408f0691795defb075912b10c9fffabf89f Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 27 Jul 2026 09:07:00 +0200 Subject: [PATCH] Address code-review findings; make detection limits detector-driven One changeset, developed together in response to a review of this branch, so the files carry several of the changes at once. Full test suite passes (733 cases). Spot finding - Split ImageSpotFinder into Detect() (flag strong pixels - the expensive per-pixel pass) and ExtractSpots() (CCL + min/max-pix + resolution mask), with Run() = both. The per-image min-pix escalation now detects ONCE and repeats only the cheap extraction, instead of re-running the whole finder four times per frame as it did on the default path. It also keeps the winning attempt's spot list rather than re-extracting it, so the frame that is integrated is exactly the frame that was scored - which a GPU re-extract could not guarantee (float atomic ordering). - spot_finding_time_s no longer swallows indexing time, and indexing_time_s now sums every escalation call instead of reporting only the last. Detection limits follow the detector - The azimuthal-integration upper q and the spot-finding high-resolution limit are now std::optional, in the C++ structs AND in the OpenAPI schema, and resolve to the detector's own maximum (DiffractionExperiment::GetDetectorMaxQ_ recipA). Adaptive detection reads a pixel's ring from the azimuthal bins, so a pixel outside that q range could never be strong - the integration range silently bounded what detection could see, regardless of the requested resolution limit. Regenerated the C++ and TypeScript clients; the viewer and the web frontend each gained a "to detector edge" switch. Detection defaults are now per workflow (measured, not assumed) - Stills: adaptive detection, min-pix chosen per image, no resolution clipping. - Rotation: fixed-threshold finder, min-pix 2, 1.5 A limit. On a 33-crystal rotation battery, adaptive detection helped four hard crystals but deterministically broke three (a lost space group, a halved indexing rate, a collapsed merge), and the detector-edge limit cost indexing on a strong rotation set (100.0 -> 96.8%). Each is still overridable by its flag, and --no-adaptive-spots is new. Indexer seed escalation - Stop escalating once a seed's lattice explains >= 90% of the seed spots. Previously any frame with >= 80 spots always paid three indexer calls, online broker included. Merge-consistency filter - --min-image-cc gated on a per-image CC computed BEFORE the stills partiality post-refinement and never refreshed; the refiner now recomputes it, so the reported CC describes the data that are actually merged. - Replaced the per-call cc_mask argument with one MergeOnTheFly flag, so the merge, the error model and MergeStats can no longer disagree about which images are in (the --scale path merged unfiltered while its statistics were filtered). Per-image B-factor refinement (-B) removed - Measured on four serial-stills datasets: it is a no-op where the per-image fit is well conditioned and actively harmful where it is not (CC1/2 -8.1, R_meas +23.2 on the weakest large-cell set, whose fits hit their [-50, 200] bounds on 14-25% of images). It had also been silently DISCARDED since the partiality post-refinement landed - reported but not applied. Rather than fix and keep a knob with no demonstrated benefit, the flag and the whole image_scale_b_factor chain are gone: setting, scaling fit, message field, CBOR, HDF5 write and read-back, per-image plot, OpenAPI enum, viewer column and checkbox, docs. ScaleOnTheFly no longer needs Ceres at all - the fit is a linear IRLS. (The Wilson per-image b_factor is a different quantity and stays.) Stills partiality width now fits both of its components - sigma^2 = gamma0^2 + (gamma_e*d*)^2 instead of a purely angular gamma_e*d* with gamma0 pinned to 0. Fitted per crystal by least squares of dist_ewald^2 on d*^2. The angular-only width is fitted over a d*^2-dense population, so it was pinned by the high-resolution edge and collapsed at low d*: median partiality 0.008 beyond 13 A for reflections that were plainly recorded, 55% of them under the merge's partiality floor, and the survivors divided by those values - which inflated the merged low-resolution intensity scale 3.6x (~ +9 A^2 of apparent B). Measured on 5000 stills: the ramp flattens to 0.89x, no observation is dropped any more (701750 -> 716811), shell-mean CC1/2 and R-free improve slightly. Note CC1/2, R_meas, completeness and a B-refining R-free are all blind to that ramp, which is why it survived earlier validation; the cost is high-resolution R_meas (98.5 -> 101.9 shell-averaged). Removed dead code from add-then-remove churn - Prediction-time "still partiality" (unreachable: no setter), the phantom IndexingSettings::min_indexed_spot_fraction knob (getter, no setter - now the constant it always was), StillsPartialityRefine's caller-less Settings constructor and its reference to a long-gone env var, ProcessImage's unread bool return, an unused include, and a dead viewer overlay hook. Also - Viewer: the magnifier compared a QImage with itself, so its scene rect was set once ever and it could not pan into a larger dataset; the hover tail timer could fire after leaveEvent and resurrect the resolution readout outside the image. - update_version.sh regenerated the frontend lock file BEFORE bumping the version (every release shipped an off-by-one lock), and did git rm/git add on a path that has not existed since the client moved to src/client - with no set -e, both failed silently. - fpga/pcie_driver/postinstall.sh tested "[ ! occurrences > 0 ]", which is a redirect, not a test, so dkms add never ran. - Unit tests for the adaptive-threshold host functions, which had none. Co-Authored-By: Claude Opus 5 (1M context) --- acquisition_device/AcquisitionDevice.cpp | 2 +- broker/OpenAPIConvert.cpp | 14 +- broker/gen/model/Azim_int_settings.cpp | 27 +++- broker/gen/model/Azim_int_settings.h | 6 +- broker/gen/model/Spot_finding_settings.cpp | 23 ++- broker/gen/model/Spot_finding_settings.h | 6 +- broker/jfjoch_api.yaml | 11 +- common/AzimuthalIntegrationSettings.cpp | 31 +++- common/AzimuthalIntegrationSettings.h | 11 +- common/BraggIntegrationSettings.cpp | 9 -- common/BraggIntegrationSettings.h | 3 - common/DiffractionExperiment.cpp | 40 ++++- common/DiffractionExperiment.h | 8 +- common/IndexingSettings.cpp | 4 - common/IndexingSettings.h | 5 - common/JFJochMessages.h | 2 - common/JFJochReceiverPlots.cpp | 8 - common/JFJochReceiverPlots.h | 1 - common/Plot.h | 2 +- common/ScalingSettings.cpp | 17 --- common/ScalingSettings.h | 9 -- docs/CBOR.md | 2 - docs/CPU_DATA_ANALYSIS.md | 16 +- docs/HDF5.md | 1 - docs/RUGNUX.md | 21 ++- fpga/pcie_driver/postinstall.sh | 2 +- frame_serialize/CBORStream2Deserializer.cpp | 4 - frame_serialize/CBORStream2Serializer.cpp | 2 - frontend/package-lock.json | 4 +- frontend/src/client/types.gen.ts | 19 ++- frontend/src/client/zod.gen.ts | 7 +- frontend/src/components/AzIntSettings.tsx | 19 ++- .../src/components/DataProcessingSettings.tsx | 21 ++- image_analysis/IndexAndRefine.cpp | 14 +- image_analysis/IndexAndRefine.h | 2 +- image_analysis/IntegrationOutcome.h | 1 - image_analysis/MXAnalysisWithoutFPGA.cpp | 55 ++++--- image_analysis/MXAnalysisWithoutFPGA.h | 4 +- .../bragg_prediction/BraggPrediction.cpp | 9 +- .../bragg_prediction/BraggPrediction.h | 7 - .../bragg_prediction/BraggPredictionGPU.cu | 17 +-- .../bragg_prediction/BraggPredictionGPU.h | 2 - image_analysis/indexing/AnalyzeIndexing.cpp | 5 +- .../rotation_indexer/RotationIndexer.cpp | 7 +- image_analysis/scale_merge/Merge.cpp | 79 ++++++++-- image_analysis/scale_merge/Merge.h | 23 ++- .../scale_merge/RotationScaleMerge.cpp | 3 +- image_analysis/scale_merge/ScaleOnTheFly.cpp | 140 ++---------------- image_analysis/scale_merge/ScaleOnTheFly.h | 7 +- image_analysis/scale_merge/ScalingResult.cpp | 6 +- image_analysis/scale_merge/ScalingResult.h | 1 - .../scale_merge/StillsPartialityRefine.cpp | 84 +++++++---- .../scale_merge/StillsPartialityRefine.h | 6 +- .../spot_finding/AdaptiveSpotFinderCPU.cpp | 16 +- .../spot_finding/AdaptiveSpotFinderCPU.h | 4 +- .../spot_finding/AdaptiveSpotFinderGPU.cu | 17 +-- .../spot_finding/AdaptiveSpotFinderGPU.h | 6 +- .../spot_finding/DetModuleSpotFinder_cpu.h | 2 +- .../spot_finding/ImageSpotFinder.cpp | 7 + image_analysis/spot_finding/ImageSpotFinder.h | 10 +- .../spot_finding/ImageSpotFinderCPU.cpp | 7 +- .../spot_finding/ImageSpotFinderCPU.h | 2 +- .../spot_finding/ImageSpotFinderGPU.cu | 4 +- .../spot_finding/ImageSpotFinderGPU.h | 2 +- .../spot_finding/SpotFindingSettings.h | 5 +- image_analysis/spot_finding/SpotUtils.cpp | 3 +- reader/HDF5MetadataSource.cpp | 3 - reader/JFJochReaderDataset.h | 1 - rugnux/Rugnux.cpp | 29 ++-- rugnux/RugnuxCommandLine.cpp | 14 +- rugnux/rugnux_cli.cpp | 72 +++++---- tests/AdaptiveThresholdTest.cpp | 100 +++++++++++++ tests/CMakeLists.txt | 2 + tests/JFJochReaderTest.cpp | 6 + tests/XDSPluginTest.cpp | 2 + update_version.sh | 18 ++- viewer/JFJochHttpReader.cpp | 1 - viewer/JFJochProcessController.cpp | 1 - viewer/JFJochViewerWindow.cpp | 7 +- .../image_viewer/JFJochDiffractionImage.cpp | 5 - viewer/image_viewer/JFJochDiffractionImage.h | 1 - viewer/image_viewer/JFJochFollowerImage.cpp | 8 +- viewer/image_viewer/JFJochFollowerImage.h | 1 + viewer/image_viewer/JFJochImage.cpp | 10 +- viewer/image_viewer/JFJochImage.h | 2 +- viewer/widgets/JFJochViewerSettingsDock.cpp | 82 +++++++--- viewer/widgets/JFJochViewerSettingsDock.h | 6 + .../windows/JFJochViewerImageListWindow.cpp | 23 +-- viewer/windows/JFJochViewerImageListWindow.h | 3 +- writer/HDF5DataFilePluginMX.cpp | 3 - writer/HDF5DataFilePluginMX.h | 1 - writer/HDF5NXmx.cpp | 1 - 92 files changed, 762 insertions(+), 554 deletions(-) create mode 100644 tests/AdaptiveThresholdTest.cpp diff --git a/acquisition_device/AcquisitionDevice.cpp b/acquisition_device/AcquisitionDevice.cpp index ed410ab9..00fecac8 100644 --- a/acquisition_device/AcquisitionDevice.cpp +++ b/acquisition_device/AcquisitionDevice.cpp @@ -323,7 +323,7 @@ void AcquisitionDevice::SetSpotFinderParameters(const SpotFindingSettings &setti fpga_parameters.snr_threshold = settings.signal_to_noise_threshold; 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_d = settings.high_resolution_limit.value_or(0.0f); 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 4783dcb2..3d6378e0 100644 --- a/broker/OpenAPIConvert.cpp +++ b/broker/OpenAPIConvert.cpp @@ -16,7 +16,8 @@ SpotFindingSettings Convert(const org::openapitools::server::model::Spot_finding ret.photon_count_threshold = input.getPhotonCountThreshold(); ret.min_pix_per_spot = input.getMinPixPerSpot(); ret.max_pix_per_spot = input.getMaxPixPerSpot(); - ret.high_resolution_limit = input.getHighResolutionLimit(); + if (input.highResolutionLimitIsSet()) + ret.high_resolution_limit = input.getHighResolutionLimit(); ret.low_resolution_limit = input.getLowResolutionLimit(); ret.enable = input.isEnable(); ret.indexing = input.isIndexing(); @@ -34,7 +35,8 @@ org::openapitools::server::model::Spot_finding_settings Convert(const SpotFindin ret.setPhotonCountThreshold(input.photon_count_threshold); ret.setMinPixPerSpot(input.min_pix_per_spot.value_or(2)); ret.setMaxPixPerSpot(input.max_pix_per_spot); - ret.setHighResolutionLimit(input.high_resolution_limit); + if (input.high_resolution_limit.has_value()) + ret.setHighResolutionLimit(input.high_resolution_limit.value()); ret.setLowResolutionLimit(input.low_resolution_limit); ret.setEnable(input.enable); ret.setIndexing(input.indexing); @@ -437,7 +439,9 @@ AzimuthalIntegrationSettings Convert(const org::openapitools::server::model::Azi ret.SolidAngleCorrection(input.isSolidAngleCorr()); ret.PolarizationCorrection(input.isPolarizationCorr()); ret.QSpacing_recipA(input.getQSpacing()); - ret.QRange_recipA(input.getLowQRecipA(), input.getHighQRecipA()); + ret.QRange_recipA(input.getLowQRecipA(), + input.highQRecipAIsSet() ? std::optional(input.getHighQRecipA()) + : std::nullopt); ret.AzimuthalBinCount(input.getAzimuthalBins()); ret.ForceCPUinFPGAWorkflow(input.isForceCpu()); return ret; @@ -447,7 +451,8 @@ org::openapitools::server::model::Azim_int_settings Convert(const AzimuthalInteg org::openapitools::server::model::Azim_int_settings ret{}; ret.setSolidAngleCorr(settings.IsSolidAngleCorrection()); ret.setPolarizationCorr(settings.IsPolarizationCorrection()); - ret.setHighQRecipA(settings.GetHighQ_recipA()); + if (const auto high_q = settings.GetRequestedHighQ_recipA()) + ret.setHighQRecipA(high_q.value()); ret.setLowQRecipA(settings.GetLowQ_recipA()); ret.setQSpacing(settings.GetQSpacing_recipA()); ret.setAzimuthalBins(settings.GetAzimuthalBinCount()); @@ -910,7 +915,6 @@ PlotType ConvertPlotType(const std::optional& input) { if (input == "integrated_reflections") return PlotType::IntegratedReflections; if (input == "image_scale_factor") return PlotType::ImageScaleFactor; if (input == "image_scale_cc") return PlotType::ImageScaleCC; - if (input == "image_scale_b") return PlotType::ImageScaleBFactor; if (input == "compression_ratio") return PlotType::CompressionRatio; if (input == "indexing_lattice_count") return PlotType::IndexingLatticeCount; throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, diff --git a/broker/gen/model/Azim_int_settings.cpp b/broker/gen/model/Azim_int_settings.cpp index f52d9820..1386b470 100644 --- a/broker/gen/model/Azim_int_settings.cpp +++ b/broker/gen/model/Azim_int_settings.cpp @@ -24,6 +24,7 @@ Azim_int_settings::Azim_int_settings() m_Polarization_corr = true; m_Solid_angle_corr = true; m_High_q_recipA = 0.0f; + m_High_q_recipAIsSet = false; m_Low_q_recipA = 0.0f; m_Q_spacing = 0.0f; m_Azimuthal_bins = 1L; @@ -53,8 +54,8 @@ bool Azim_int_settings::validate(std::stringstream& msg, const std::string& path const std::string _pathPrefix = pathPrefix.empty() ? "Azim_int_settings" : pathPrefix; - - /* High_q_recipA */ { + if (highQRecipAIsSet()) + { const float& value = m_High_q_recipA; const std::string currentValuePath = _pathPrefix + ".highQRecipA"; @@ -138,8 +139,8 @@ bool Azim_int_settings::operator==(const Azim_int_settings& rhs) const (isSolidAngleCorr() == rhs.isSolidAngleCorr()) && - (getHighQRecipA() == rhs.getHighQRecipA()) - && + + ((!highQRecipAIsSet() && !rhs.highQRecipAIsSet()) || (highQRecipAIsSet() && rhs.highQRecipAIsSet() && getHighQRecipA() == rhs.getHighQRecipA())) && (getLowQRecipA() == rhs.getLowQRecipA()) && @@ -166,7 +167,8 @@ void to_json(nlohmann::json& j, const Azim_int_settings& o) j = nlohmann::json::object(); j["polarization_corr"] = o.m_Polarization_corr; j["solid_angle_corr"] = o.m_Solid_angle_corr; - j["high_q_recipA"] = o.m_High_q_recipA; + if(o.highQRecipAIsSet()) + j["high_q_recipA"] = o.m_High_q_recipA; j["low_q_recipA"] = o.m_Low_q_recipA; j["q_spacing"] = o.m_Q_spacing; if(o.azimuthalBinsIsSet()) @@ -180,7 +182,11 @@ void from_json(const nlohmann::json& j, Azim_int_settings& o) { j.at("polarization_corr").get_to(o.m_Polarization_corr); j.at("solid_angle_corr").get_to(o.m_Solid_angle_corr); - j.at("high_q_recipA").get_to(o.m_High_q_recipA); + if(j.find("high_q_recipA") != j.end()) + { + j.at("high_q_recipA").get_to(o.m_High_q_recipA); + o.m_High_q_recipAIsSet = true; + } j.at("low_q_recipA").get_to(o.m_Low_q_recipA); j.at("q_spacing").get_to(o.m_Q_spacing); if(j.find("azimuthal_bins") != j.end()) @@ -219,6 +225,15 @@ float Azim_int_settings::getHighQRecipA() const void Azim_int_settings::setHighQRecipA(float const value) { m_High_q_recipA = value; + m_High_q_recipAIsSet = true; +} +bool Azim_int_settings::highQRecipAIsSet() const +{ + return m_High_q_recipAIsSet; +} +void Azim_int_settings::unsetHigh_q_recipA() +{ + m_High_q_recipAIsSet = false; } float Azim_int_settings::getLowQRecipA() const { diff --git a/broker/gen/model/Azim_int_settings.h b/broker/gen/model/Azim_int_settings.h index a412f93f..51889fd7 100644 --- a/broker/gen/model/Azim_int_settings.h +++ b/broker/gen/model/Azim_int_settings.h @@ -68,10 +68,12 @@ public: bool isSolidAngleCorr() const; void setSolidAngleCorr(bool const value); /// - /// + /// Upper q limit of the azimuthal integration [1/Angstrom]. Optional: if omitted, the integration (and the adaptive spot detection that shares these q bins) extends to the highest q the detector reaches. /// float getHighQRecipA() const; void setHighQRecipA(float const value); + bool highQRecipAIsSet() const; + void unsetHigh_q_recipA(); /// /// /// @@ -105,7 +107,7 @@ protected: bool m_Solid_angle_corr; float m_High_q_recipA; - + bool m_High_q_recipAIsSet; float m_Low_q_recipA; float m_Q_spacing; diff --git a/broker/gen/model/Spot_finding_settings.cpp b/broker/gen/model/Spot_finding_settings.cpp index 7f178aca..b6cc66b7 100644 --- a/broker/gen/model/Spot_finding_settings.cpp +++ b/broker/gen/model/Spot_finding_settings.cpp @@ -28,6 +28,7 @@ Spot_finding_settings::Spot_finding_settings() m_Min_pix_per_spot = 0L; m_Max_pix_per_spot = 0L; m_High_resolution_limit = 0.0f; + m_High_resolution_limitIsSet = false; m_Low_resolution_limit = 0.0f; m_High_resolution_limit_for_spot_count_low_res = 0.0f; m_Quick_integration = false; @@ -196,8 +197,8 @@ bool Spot_finding_settings::operator==(const Spot_finding_settings& rhs) const (getMaxPixPerSpot() == rhs.getMaxPixPerSpot()) && - (getHighResolutionLimit() == rhs.getHighResolutionLimit()) - && + + ((!highResolutionLimitIsSet() && !rhs.highResolutionLimitIsSet()) || (highResolutionLimitIsSet() && rhs.highResolutionLimitIsSet() && getHighResolutionLimit() == rhs.getHighResolutionLimit())) && (getLowResolutionLimit() == rhs.getLowResolutionLimit()) && @@ -231,7 +232,8 @@ void to_json(nlohmann::json& j, const Spot_finding_settings& o) j["photon_count_threshold"] = o.m_Photon_count_threshold; j["min_pix_per_spot"] = o.m_Min_pix_per_spot; j["max_pix_per_spot"] = o.m_Max_pix_per_spot; - j["high_resolution_limit"] = o.m_High_resolution_limit; + if(o.highResolutionLimitIsSet()) + j["high_resolution_limit"] = o.m_High_resolution_limit; j["low_resolution_limit"] = o.m_Low_resolution_limit; j["high_resolution_limit_for_spot_count_low_res"] = o.m_High_resolution_limit_for_spot_count_low_res; j["quick_integration"] = o.m_Quick_integration; @@ -249,7 +251,11 @@ void from_json(const nlohmann::json& j, Spot_finding_settings& o) j.at("photon_count_threshold").get_to(o.m_Photon_count_threshold); j.at("min_pix_per_spot").get_to(o.m_Min_pix_per_spot); j.at("max_pix_per_spot").get_to(o.m_Max_pix_per_spot); - j.at("high_resolution_limit").get_to(o.m_High_resolution_limit); + if(j.find("high_resolution_limit") != j.end()) + { + j.at("high_resolution_limit").get_to(o.m_High_resolution_limit); + o.m_High_resolution_limitIsSet = true; + } j.at("low_resolution_limit").get_to(o.m_Low_resolution_limit); j.at("high_resolution_limit_for_spot_count_low_res").get_to(o.m_High_resolution_limit_for_spot_count_low_res); j.at("quick_integration").get_to(o.m_Quick_integration); @@ -317,6 +323,15 @@ float Spot_finding_settings::getHighResolutionLimit() const void Spot_finding_settings::setHighResolutionLimit(float const value) { m_High_resolution_limit = value; + m_High_resolution_limitIsSet = true; +} +bool Spot_finding_settings::highResolutionLimitIsSet() const +{ + return m_High_resolution_limitIsSet; +} +void Spot_finding_settings::unsetHigh_resolution_limit() +{ + m_High_resolution_limitIsSet = false; } float Spot_finding_settings::getLowResolutionLimit() const { diff --git a/broker/gen/model/Spot_finding_settings.h b/broker/gen/model/Spot_finding_settings.h index 4eb5f17c..6acf3381 100644 --- a/broker/gen/model/Spot_finding_settings.h +++ b/broker/gen/model/Spot_finding_settings.h @@ -88,10 +88,12 @@ public: int64_t getMaxPixPerSpot() const; void setMaxPixPerSpot(int64_t const value); /// - /// High resolution limit for spot finding [Angstrom] + /// High resolution limit for spot finding [Angstrom]. Optional: if omitted, spot finding extends as far as the detector reaches, i.e. the detection is not clipped in resolution. /// float getHighResolutionLimit() const; void setHighResolutionLimit(float const value); + bool highResolutionLimitIsSet() const; + void unsetHigh_resolution_limit(); /// /// Low resolution limit for spot finding [Angstrom] /// @@ -136,7 +138,7 @@ protected: int64_t m_Max_pix_per_spot; float m_High_resolution_limit; - + bool m_High_resolution_limitIsSet; float m_Low_resolution_limit; float m_High_resolution_limit_for_spot_count_low_res; diff --git a/broker/jfjoch_api.yaml b/broker/jfjoch_api.yaml index 67a4e5fc..4385fe5a 100644 --- a/broker/jfjoch_api.yaml +++ b/broker/jfjoch_api.yaml @@ -119,7 +119,6 @@ components: - integrated_reflections - image_scale_factor - image_scale_cc - - image_scale_b - compression_ratio - ice_ring_score roi: @@ -1039,7 +1038,6 @@ components: - photon_count_threshold - max_pix_per_spot - min_pix_per_spot - - high_resolution_limit - low_resolution_limit - quick_integration - high_resolution_limit_for_spot_count_low_res @@ -1075,7 +1073,9 @@ components: high_resolution_limit: type: number format: float - description: High resolution limit for spot finding [Angstrom] + description: | + High resolution limit for spot finding [Angstrom]. Optional: if omitted, spot finding extends + as far as the detector reaches, i.e. the detection is not clipped in resolution. low_resolution_limit: type: number format: float @@ -1115,7 +1115,6 @@ components: required: - solid_angle_corr - polarization_corr - - high_q_recipA - low_q_recipA - q_spacing properties: @@ -1132,6 +1131,10 @@ components: minimum: 2e-5 maximum: 10.0 format: float + description: | + Upper q limit of the azimuthal integration [1/Angstrom]. Optional: if omitted, the integration + (and the adaptive spot detection that shares these q bins) extends to the highest q the + detector reaches. low_q_recipA: type: number format: float diff --git a/common/AzimuthalIntegrationSettings.cpp b/common/AzimuthalIntegrationSettings.cpp index 76ccf4ff..de5504c9 100644 --- a/common/AzimuthalIntegrationSettings.cpp +++ b/common/AzimuthalIntegrationSettings.cpp @@ -1,6 +1,7 @@ // SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute // SPDX-License-Identifier: GPL-3.0-only +#include #include #include "AzimuthalIntegrationSettings.h" @@ -19,21 +20,33 @@ AzimuthalIntegrationSettings &AzimuthalIntegrationSettings::SolidAngleCorrection return *this; } -AzimuthalIntegrationSettings &AzimuthalIntegrationSettings::QRange_recipA(float low, float high) { +AzimuthalIntegrationSettings &AzimuthalIntegrationSettings::QRange_recipA(float low, std::optional high) { check_finite("Low Q for azimuthal integration", low); - check_finite("High Q for azimuthal integration", high); - check_max("High Q for azimuthal integration", high, maxQ_recipA); check_min("Low Q for azimuthal integration", low, minQ_recipA); - if (high <= low) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - "High Q must be higher than low Q"); + if (high.has_value()) { + check_finite("High Q for azimuthal integration", *high); + check_max("High Q for azimuthal integration", *high, maxQ_recipA); + if (*high <= low) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "High Q must be higher than low Q"); + } - high_q_recipA = high; + requested_high_q_recipA = high; low_q_recipA = low; + // Until ResolveHighQ runs, an unset limit keeps the value the bins were last built from. + high_q_recipA = high.value_or(high_q_recipA); UpdateBinCount(); return *this; } +void AzimuthalIntegrationSettings::ResolveHighQ(float detector_max_q_recipA) { + if (requested_high_q_recipA.has_value()) + return; + + high_q_recipA = std::clamp(detector_max_q_recipA, low_q_recipA + q_spacing, maxQ_recipA); + UpdateBinCount(); +} + AzimuthalIntegrationSettings &AzimuthalIntegrationSettings::QSpacing_recipA(float input) { check_finite("Q spacing for azimuthal integration", input); check_min("Q spacing for azimuthal integration", input, minQ_recipA); @@ -50,6 +63,10 @@ float AzimuthalIntegrationSettings::GetHighQ_recipA() const { return high_q_recipA; } +std::optional AzimuthalIntegrationSettings::GetRequestedHighQ_recipA() const { + return requested_high_q_recipA; +} + float AzimuthalIntegrationSettings::GetLowQ_recipA() const { return low_q_recipA; } diff --git a/common/AzimuthalIntegrationSettings.h b/common/AzimuthalIntegrationSettings.h index cc750657..63835231 100644 --- a/common/AzimuthalIntegrationSettings.h +++ b/common/AzimuthalIntegrationSettings.h @@ -13,6 +13,12 @@ class AzimuthalIntegrationSettings { bool solid_angle_correction = true; bool polarization_correction = true; + // Requested upper q limit. Unset means "as far as the detector reaches": DiffractionExperiment + // resolves it from the geometry (ResolveHighQ) whenever it hands these settings out, so + // high_q_recipA below - what the bins are built from - is always a concrete number. Not clipping + // detection at an arbitrary default matters for the adaptive spot finder, which bins pixels through + // this same q range and cannot see a pixel that falls outside it. + std::optional requested_high_q_recipA; float high_q_recipA = 5.0; float low_q_recipA = 0.1; float bkg_estimate_high_q_recipA = 2.0f * PI / 3.0; @@ -31,7 +37,9 @@ public: AzimuthalIntegrationSettings(); AzimuthalIntegrationSettings& SolidAngleCorrection(bool input); AzimuthalIntegrationSettings& PolarizationCorrection(bool input); - AzimuthalIntegrationSettings& QRange_recipA(float low, float high); + AzimuthalIntegrationSettings& QRange_recipA(float low, std::optional high); + // Substitute the detector's own maximum q for an unset high q. No-op if one was requested. + void ResolveHighQ(float detector_max_q_recipA); AzimuthalIntegrationSettings& QSpacing_recipA(float input); AzimuthalIntegrationSettings& BkgEstimateQRange_recipA(float low, float high); AzimuthalIntegrationSettings& AzimuthalBinCount(int32_t input); @@ -40,6 +48,7 @@ public: [[nodiscard]] bool IsSolidAngleCorrection() const; [[nodiscard]] bool IsPolarizationCorrection() const; [[nodiscard]] float GetHighQ_recipA() const; + [[nodiscard]] std::optional GetRequestedHighQ_recipA() const; [[nodiscard]] float GetLowQ_recipA() const; [[nodiscard]] float GetQSpacing_recipA() const; diff --git a/common/BraggIntegrationSettings.cpp b/common/BraggIntegrationSettings.cpp index 338743d8..aa8ae524 100644 --- a/common/BraggIntegrationSettings.cpp +++ b/common/BraggIntegrationSettings.cpp @@ -97,15 +97,6 @@ float BraggIntegrationSettings::GetMinimumSigmaInRegardsToI() const { return minimum_sigma_in_regards_to_i; } -BraggIntegrationSettings &BraggIntegrationSettings::StillPartiality(bool input) { - still_partiality = input; - return *this; -} - -bool BraggIntegrationSettings::GetStillPartiality() const { - return still_partiality; -} - BraggIntegrationSettings &BraggIntegrationSettings::BackgroundTrimFraction(float input) { check_finite("Background trim fraction", input); check_min("Background trim fraction", input, 0.0); diff --git a/common/BraggIntegrationSettings.h b/common/BraggIntegrationSettings.h index 36b31b9f..974dc47c 100644 --- a/common/BraggIntegrationSettings.h +++ b/common/BraggIntegrationSettings.h @@ -20,7 +20,6 @@ class BraggIntegrationSettings { float d_min_limit_A = 1.0; std::optional fixed_profile_radius; float minimum_sigma_in_regards_to_i = 0.02; - bool still_partiality = false; // experimental stills excitation-error partiality (rugnux --still-partiality) // Symmetric trimmed-mean fraction for the r2..r3 background ring: drop the lowest and highest this // fraction of ring pixels before averaging. Resists the high-side contamination (neighbour-spot // wings, tails, zingers) that biases the plain ring mean up and makes it over-subtract weak @@ -36,7 +35,6 @@ public: BraggIntegrationSettings& DMinLimit_A(float input); BraggIntegrationSettings& FixedProfileRadius_recipA(std::optional input); BraggIntegrationSettings& Integrator(IntegratorMode input); - BraggIntegrationSettings& StillPartiality(bool input); BraggIntegrationSettings& BackgroundTrimFraction(float input); @@ -48,6 +46,5 @@ public: [[nodiscard]] float GetDMinLimit_A() const; [[nodiscard]] float GetMinimumSigmaInRegardsToI() const; - [[nodiscard]] bool GetStillPartiality() const; [[nodiscard]] float GetBackgroundTrimFraction() const; }; diff --git a/common/DiffractionExperiment.cpp b/common/DiffractionExperiment.cpp index 38f50880..21068756 100644 --- a/common/DiffractionExperiment.cpp +++ b/common/DiffractionExperiment.cpp @@ -1,6 +1,7 @@ // SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute // SPDX-License-Identifier: GPL-3.0-only +#include #include #include @@ -164,7 +165,7 @@ DiffractionExperiment &DiffractionExperiment::MaskChipEdges(bool input) { return *this; } -DiffractionExperiment &DiffractionExperiment::QRangeForAzimInt_recipA(float low, float high) { +DiffractionExperiment &DiffractionExperiment::QRangeForAzimInt_recipA(float low, std::optional high) { az_integration_settings.QRange_recipA(low, high); return *this; } @@ -581,7 +582,7 @@ float DiffractionExperiment::GetLowQForAzimInt_recipA() const { return az_integration_settings.GetLowQ_recipA(); } float DiffractionExperiment::GetHighQForAzimInt_recipA() const { - return az_integration_settings.GetHighQ_recipA(); + return GetAzimuthalIntegrationSettings().GetHighQ_recipA(); } float DiffractionExperiment::GetQSpacingForAzimInt_recipA() const { @@ -610,15 +611,16 @@ void DiffractionExperiment::CheckDataProcessingSettings(const SpotFindingSetting check_min("Photon count threshold", settings.photon_count_threshold, 0); 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); - if (settings.high_resolution_limit > 0) { - check_min("Spot finding high resolution limit", settings.high_resolution_limit, 0.5); - check_max("Spot finding high resolution limit", settings.high_resolution_limit, 50.0); + // An unset high-resolution limit means "as far as the detector reaches", so there is nothing to check. + if (settings.high_resolution_limit.value_or(0.0f) > 0) { + check_finite("Spot finding high resolution limit", *settings.high_resolution_limit); + check_min("Spot finding high resolution limit", *settings.high_resolution_limit, 0.5); + check_max("Spot finding high resolution limit", *settings.high_resolution_limit, 50.0); if (settings.low_resolution_limit > 0) { check_min("Spot finding low resolution limit", settings.low_resolution_limit, - settings.high_resolution_limit); + *settings.high_resolution_limit); } } else if (settings.low_resolution_limit > 0) { check_min("Spot finding low resolution limit", settings.low_resolution_limit, 1.0); @@ -1316,7 +1318,29 @@ DiffractionExperiment &DiffractionExperiment::ImportAzimuthalIntegrationSettings } AzimuthalIntegrationSettings DiffractionExperiment::GetAzimuthalIntegrationSettings() const { - return az_integration_settings; + // An unset high q means "as far as the detector reaches", so resolve it here, where the geometry is + // known. Everyone reads the settings through this getter, so nobody sees an unresolved q range. + AzimuthalIntegrationSettings ret = az_integration_settings; + ret.ResolveHighQ(GetDetectorMaxQ_recipA()); + return ret; +} + +float DiffractionExperiment::GetDetectorMaxResolution_A() const { + const float q = GetDetectorMaxQ_recipA(); + return q > 0.0f ? 2.0f * static_cast(PI) / q : 0.0f; +} + +float DiffractionExperiment::GetDetectorMaxQ_recipA() const { + const DiffractionGeometry geom = GetDiffractionGeometry(); + const auto width = static_cast(GetXPixelsNumConv()); + const auto height = static_cast(GetYPixelsNumConv()); + + // The largest scattering angle sits at one of the detector corners, wherever the beam centre is. + float q = 0.0f; + for (const float x: {0.0f, width}) + for (const float y: {0.0f, height}) + q = std::max(q, geom.PxlToQ(x, y)); + return q; } DiffractionExperiment &DiffractionExperiment::PolarizationFactor(const std::optional &input) { diff --git a/common/DiffractionExperiment.h b/common/DiffractionExperiment.h index 98167523..f428961d 100644 --- a/common/DiffractionExperiment.h +++ b/common/DiffractionExperiment.h @@ -105,7 +105,7 @@ public: DiffractionExperiment& MaskModuleEdges(bool input); DiffractionExperiment& MaskChipEdges(bool input); - DiffractionExperiment& QRangeForAzimInt_recipA(float low, float high); + DiffractionExperiment& QRangeForAzimInt_recipA(float low, std::optional high); DiffractionExperiment& BkgEstimateQRange_recipA(float low, float high); DiffractionExperiment& QSpacingForAzimInt_recipA(float input); @@ -175,6 +175,12 @@ public: DiffractionExperiment& ImportAzimuthalIntegrationSettings(const AzimuthalIntegrationSettings& input); AzimuthalIntegrationSettings GetAzimuthalIntegrationSettings() const; + // Highest q (2*pi/d) any pixel of the detector reaches, from the current geometry, and the same + // limit as a resolution in Angstrom. This is what an unset azimuthal-integration high q and an unset + // spot-finding high-resolution limit resolve to. + [[nodiscard]] float GetDetectorMaxQ_recipA() const; + [[nodiscard]] float GetDetectorMaxResolution_A() const; + DiffractionExperiment& ImportBraggIntegrationSettings(const BraggIntegrationSettings& input); BraggIntegrationSettings GetBraggIntegrationSettings() const; diff --git a/common/IndexingSettings.cpp b/common/IndexingSettings.cpp index 122bd0e8..2de46bcb 100644 --- a/common/IndexingSettings.cpp +++ b/common/IndexingSettings.cpp @@ -28,10 +28,6 @@ int64_t IndexingSettings::GetViableCellMinSpots() const { return viable_cell_min_spots; } -float IndexingSettings::GetMinIndexedSpotFraction() const { - return min_indexed_spot_fraction; -} - IndexingSettings &IndexingSettings::Algorithm(IndexingAlgorithmEnum input) { switch (input) { case IndexingAlgorithmEnum::Auto: diff --git a/common/IndexingSettings.h b/common/IndexingSettings.h index 7c8ee11d..258e15ba 100644 --- a/common/IndexingSettings.h +++ b/common/IndexingSettings.h @@ -24,10 +24,6 @@ class IndexingSettings { static constexpr float unit_cell_angle_tolerance_deg = 5.0; // degree int64_t indexing_threads = 4; int64_t viable_cell_min_spots = 9; - // Minimum fraction of the in-resolution spots a candidate lattice must index to be accepted. - // Lowering it admits weaker/sparser crystals (more real ones on flooded XFEL frames, but also - // more spurious lattices that a downstream merge-consistency gate must remove). - float min_indexed_spot_fraction = 0.20f; int64_t max_extra_lattices = 2; @@ -62,7 +58,6 @@ public: IndexingSettings& MaxExtraLattices(int64_t input); [[nodiscard]] int64_t GetViableCellMinSpots() const; - [[nodiscard]] float GetMinIndexedSpotFraction() const; [[nodiscard]] IndexingAlgorithmEnum GetAlgorithm() const; [[nodiscard]] GeomRefinementAlgorithmEnum GetGeomRefinementAlgorithm() const; [[nodiscard]] float GetFFT_MaxUnitCell_A() const; diff --git a/common/JFJochMessages.h b/common/JFJochMessages.h index a884d4c9..940bd37c 100644 --- a/common/JFJochMessages.h +++ b/common/JFJochMessages.h @@ -185,7 +185,6 @@ struct DataMessage { std::optional image_scale_factor; std::optional image_scale_cc; std::optional image_scale_mosaicity; - std::optional image_scale_b_factor; std::optional compression_ratio; }; @@ -386,7 +385,6 @@ struct EndMessage { std::vector image_scale_factor; std::vector image_scale_cc; - std::vector image_scale_b_factor; std::vector image_scale_mosaicity; std::vector ice_ring_score; }; diff --git a/common/JFJochReceiverPlots.cpp b/common/JFJochReceiverPlots.cpp index 16a22b50..7802dadc 100644 --- a/common/JFJochReceiverPlots.cpp +++ b/common/JFJochReceiverPlots.cpp @@ -109,7 +109,6 @@ void JFJochReceiverPlots::Setup(const DiffractionExperiment &experiment, const A integrated_reflections.Clear(r); image_scale_factor.Clear(r); image_scale_cc.Clear(r); - image_scale_b.Clear(r); refinement_time.Clear(r); spot_finding_time.Clear(r); @@ -180,7 +179,6 @@ void JFJochReceiverPlots::Add(const DataMessage &msg, const AzimuthalIntegration indexing_solution.AddElement(msg.number, msg.indexing_result); image_scale_factor.AddElement(msg.number, msg.image_scale_factor); image_scale_cc.AddElement(msg.number, msg.image_scale_cc); - image_scale_b.AddElement(msg.number, msg.image_scale_b_factor); { std::unique_lock ul(m); @@ -345,9 +343,6 @@ MultiLinePlot JFJochReceiverPlots::GetPlots(const PlotRequest &request) { case PlotType::ImageScaleCC: ret = image_scale_cc.GetMeanPlot(nbins, start, incr, request.fill_value); break; - case PlotType::ImageScaleBFactor: - ret = image_scale_b.GetMeanPlot(nbins, start, incr, request.fill_value); - break; case PlotType::ImageScaleFactor: ret = image_scale_factor.GetMeanPlot(nbins, start, incr, request.fill_value); break; @@ -575,9 +570,6 @@ void JFJochReceiverPlots::GetPlotRaw(std::vector &v, PlotType type, const case PlotType::ImageScaleFactor: v = image_scale_factor.ExportArray(); break; - case PlotType::ImageScaleBFactor: - v = image_scale_b.ExportArray(); - break; case PlotType::CompressionRatio: v = compression_ratio.ExportArray(); break; diff --git a/common/JFJochReceiverPlots.h b/common/JFJochReceiverPlots.h index 01b7c17c..4af96d24 100644 --- a/common/JFJochReceiverPlots.h +++ b/common/JFJochReceiverPlots.h @@ -70,7 +70,6 @@ class JFJochReceiverPlots { StatusVector integrated_reflections; StatusVector image_scale_factor; StatusVector image_scale_cc; - StatusVector image_scale_b; StatusVector compression_ratio; // StatusVector objects are fully thread-safe (protected by internal mutex) diff --git a/common/Plot.h b/common/Plot.h index 4ecbe1c1..b629c1b4 100644 --- a/common/Plot.h +++ b/common/Plot.h @@ -15,7 +15,7 @@ enum class PlotType { ROISum, ROIMean, ROIMaxCount, ROIPixels, ROIWeightedX, ROIWeightedY, PacketsReceived, MaxValue, ResolutionEstimate, ProfileRadius, Mosaicity, BFactor, PixelSum, StrongPixels, RefinementBeamX, RefinementBeamY, ImageProcessingTime, IntegratedReflections, - ImageScaleFactor, ImageScaleCC, ImageScaleBFactor, CompressionRatio, IndexingLatticeCount, IceRingScore + ImageScaleFactor, ImageScaleCC, CompressionRatio, IndexingLatticeCount, IceRingScore }; enum class PlotAzintUnit { diff --git a/common/ScalingSettings.cpp b/common/ScalingSettings.cpp index fe7a8952..9a43f393 100644 --- a/common/ScalingSettings.cpp +++ b/common/ScalingSettings.cpp @@ -5,11 +5,6 @@ #include "ScalingSettings.h" -ScalingSettings& ScalingSettings::RefineB(bool input) { - refine_b = input; - return *this; -} - ScalingSettings& ScalingSettings::MergeFriedel(bool input) { merge_friedel = input; return *this; @@ -29,10 +24,6 @@ ScalingSettings& ScalingSettings::HighResolutionLimit_A(std::optional li return *this; } -bool ScalingSettings::GetRefineB() const { - return refine_b; -} - bool ScalingSettings::GetMergeFriedel() const { return merge_friedel; } @@ -50,14 +41,6 @@ std::optional ScalingSettings::GetHighResolutionLimit_A() const { return high_resolution_limit_A; } -double ScalingSettings::GetMinB() const { - return min_b; -} - -double ScalingSettings::GetMaxB() const { - return max_b; -} - double ScalingSettings::GetMinMosaicity() const { return 0.001; } diff --git a/common/ScalingSettings.h b/common/ScalingSettings.h index b372e00b..abef9dfa 100644 --- a/common/ScalingSettings.h +++ b/common/ScalingSettings.h @@ -12,10 +12,6 @@ enum class ResolutionCutoffMethod { Off, CCHalfLogistic }; class ScalingSettings { - bool refine_b = false; - double max_b = 200.0; - double min_b = -50.0; - bool refine_wedge = false; bool merge_friedel = true; @@ -86,7 +82,6 @@ class ScalingSettings { bool scaling_regularize = false; public: - ScalingSettings& RefineB(bool input); ScalingSettings& RefineRotationWedge(bool input); ScalingSettings& RotationWedgeForScaling(std::optional input); ScalingSettings& MergeFriedel(bool input); @@ -113,12 +108,8 @@ public: ScalingSettings& ResolutionCCTarget(double input); ScalingSettings& ReportShellCount(int input); - [[nodiscard]] bool GetRefineB() const; [[nodiscard]] bool GetRefineWedge() const; - [[nodiscard]] double GetMinB() const; - [[nodiscard]] double GetMaxB() const; - [[nodiscard]] double GetMinMosaicity() const; [[nodiscard]] double GetDefaultMosaicity() const; [[nodiscard]] double GetMaxMosaicity() const; diff --git a/docs/CBOR.md b/docs/CBOR.md index 5f3fd2c5..0dbcd207 100644 --- a/docs/CBOR.md +++ b/docs/CBOR.md @@ -220,7 +220,6 @@ See [DECTRIS documentation](https://github.com/dectris/documentation/tree/main/s | beam_corr_y | float | Beam center correction Y applied during processing \[pixel\] | | X | | image_scale_factor | float | Scaling result: Image scale factor (g) | | X | | image_scale_mosaicity | float | Scaling result: Image scale mosaicity \[deg\] | | X | -| image_scale_b_factor | float | Scaling result: Image scale B factor \[Angstrom^2\] | | X | | image_scale_cc | float | Scaling result: Image scale CC | | X | | adu_histogram | Array(uint64) | ADU histogram | | | | roi_integrals | object | Results of ROI calculation | | X | @@ -321,7 +320,6 @@ See [DECTRIS documentation](https://github.com/dectris/documentation/tree/main/s | niggli_class | Array(uint8) | Per-image Niggli class identifier for indexed images; 0 if unavailable | | | pixel_sum | Array(int64) | Per-image sum of all valid pixels, excluding error/saturated pixels | | | image_scale_mosaicity | Array(float) | Scaling result: Image scale mosaicity \[deg\] | | -| image_scale_b_factor | Array(float) | Scaling result: Image scale B factor \[Angstrom^2\] | | | image_scale_cc | Array(float) | Scaling result: Image scale CC | | End-message vector fields are optional. When present, they provide master-file summary data so readers can inspect scan-level and per-image analysis results without opening every linked data file. Missing optional per-image values are encoded by the producer as zero unless otherwise noted. diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index c284ac71..234dfff2 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -187,7 +187,7 @@ Special cases: ### 3.2 Adaptive (self-calibrating) detection -The local-statistics test above still needs a fixed photon/count threshold, and the right value depends on the background level, which varies between datasets — so it has to be tuned per dataset. An optional **adaptive** mode (`--adaptive-spots`) removes that tuning by deriving the threshold from each image's own noise, per resolution ring. +The local-statistics test above still needs a fixed photon/count threshold, and the right value depends on the background level, which varies between datasets — so it has to be tuned per dataset. The **adaptive** mode (`--adaptive-spots`; the default for stills in `rugnux` and in the viewer, `--no-adaptive-spots` reverts) removes that tuning by deriving the threshold from each image's own noise, per resolution ring. Rotation data keeps the fixed-threshold finder by default: across a 33-crystal rotation battery adaptive detection helped four hard crystals but deterministically broke three (a lost space group, a halved indexing rate, a collapsed merge). Pixels are binned into the same resolution rings as the azimuthal integrator (§2). For each ring a robust background is estimated in three passes: one plain pass over all valid pixels, then two $\sigma$-clipping passes that keep only pixels within $\pm 3\sigma$ of the current ring mean (removing the Bragg peaks from the background estimate). This yields a per-ring background mean $\mu_b$ and scatter $\sigma_b$. @@ -197,6 +197,8 @@ t_b = \max\!\big(\;\mu_b + z\,\sqrt{\sigma_b^2 + \sigma_\mathrm{read}^2}\;,\;\; $ where $k_\mathrm{Poisson}(\mu_b,p)$ is the smallest count whose Poisson$(\mu_b)$ upper tail is $\le p$. The Poisson arm is correct where the background is countable (a bright low-resolution ring gets a high threshold); the Gaussian arm — floored by a detector-level excess-noise constant $\sigma_\mathrm{read}$ — takes over on near-empty high-resolution rings, where the Poisson arm degenerates to "one photon is significant" and would flood. The operating point $p = E/N$ is set from a single portable knob $E$, the expected number of false pixels tolerated per frame (`--spot-false-pixels`, default 100), with $N$ the number of valid pixels. Because $p$ and every $\mu_b,\sigma_b$ come from the image itself, the same $E$ lands a sensible photon threshold on strong and weak datasets alike, with no per-dataset tuning. Rings too sparse to characterise (fewer than ~40 pixels) fall back to a whole-frame background. A pixel is strong when $v_i \ge t_b$ for its ring (saturated pixels are still forced strong); the strong pixels then feed the same CCL stage (§3.4). The signal-to-noise and photon-count criteria of §3.1 are not used in this mode. +Because detection reads the pixel's ring, a pixel that falls outside the azimuthal-integration $q$ range has no ring and can never be strong: the integration range bounds what adaptive detection can see. Both upper limits are therefore optional and default to the detector itself — the azimuthal integration runs to the highest $q$ any pixel of the detector reaches (`--azim-max-q` unset), and for stills, spot finding is not clipped in resolution (`--spot-high-resolution` unset). Setting either one narrows detection accordingly. Rotation data keeps a 1.5 Å spot-finding limit by default: the extra high-resolution spots the detector-edge limit admits are mostly noise there and cost indexing (measured 100.0 → 96.8 % on a strong rotation set). + **Fused GPU engine.** The per-ring reduction the adaptive threshold needs is the *same* reduction the azimuthal integrator performs. On the GPU path the two are fused into a single image pass (`AdaptiveSpotFinderGPU`): one reduction accumulates the corrected per-ring sums for the azimuthal profile (§2) *and* the raw per-ring statistics for the threshold, after which a light kernel flags the strong pixels. This replaces both the separate azimuthal-integration pass and the host-side adaptive spot-finding pass with one GPU pass — on a ~4.5 MP detector it runs in ~1 ms/frame versus ~40 ms for the CPU adaptive finder, and produces an identical spot list and azimuthal profile. It is enabled by default in the offline `rugnux` path and the interactive viewer; the online receiver keeps the CPU adaptive finder. ### 3.3 Resolution and ice-ring handling @@ -225,7 +227,7 @@ The minimum-pixels-per-spot filter (§3.4) trades sensitivity against noise: a s $$ \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. +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, strong-pixel detection runs **once** per frame and the three attempts only repeat the cheap CCL and spot-level filter, not the pixel reduction; the azimuthal profile is the one that single detection pass computed. The winning attempt's spot list is kept rather than re-extracted, so the frame that is integrated is exactly the frame that was scored. 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. --- @@ -546,7 +548,7 @@ The partiality applied is fixed by the data type and scaling stage, not chosen f 2. **Unity** ($P_{ij}=1$): used for the scale-on-fulls refit (§10.6), where each observation is already a complete reflection. -3. **Fixed**: use the per-reflection partiality carried from prediction. Still/serial images are predicted with $P=1$ by default, so their scaling is effectively unity/fixed. An optional excitation-error still-partiality model (`--still-partiality`) instead weights each stills reflection by a Gaussian $\exp(-\Delta_\mathrm{Ewald}^2/2\sigma^2)$ in its distance from the Ewald sphere, with a companion merge term (`--partiality-uncertainty`) that adds an intensity- and $(1-P)$-proportional sigma to down-weight the least-complete reflections. +3. **Fixed**: use the per-reflection partiality carried from prediction. Still/serial images are predicted with $P=1$, so a single-pass stills scale is effectively unity/fixed — which is exactly what `--simple-stills` keeps. By default the stills path instead **post-refines a physical partiality**: a small per-crystal orientation tilt $(\delta\psi_x,\delta\psi_y)$ about the two axes perpendicular to the beam is refined against the running merge, and every reflection's partiality is then recomputed analytically from the refined lattice through its excitation error $\Delta_\mathrm{Ewald}=\big|\,|\mathbf{q}+\mathbf{S}_0|-1/\lambda\,\big|$ and an angular width. A tilt moves reflections on opposite sides of the Ewald sphere in opposite directions, so it reshapes the *spatial* pattern of partialities in a way the per-image scale $G$ cannot mimic — which is why it succeeds where a freely-fitted scalar partiality width simply collapses into $G$. Nothing is re-integrated (the integrated intensities are fixed); the tilt is hard-bounded at about 1° and held by a soft prior, so it stays inert on sparse or weak crystals. The cycle is merge → per-crystal tilt refinement (with $G$ profiled out by the same robust Cauchy IRLS used for the per-frame scales, §10.3) → recompute $P$ → re-merge, repeated a few times. Reflections below a minimum partiality can be rejected from merging to avoid unstable corrections. @@ -568,6 +570,8 @@ I_h = \frac{\sum_j w_j I^{\mathrm{corr}}_{ij}}{\sum_j w_j},\qquad w_j = \frac{1}{(\sigma^{\mathrm{corr}}_{ij})^2}. $ +The weights use an **expected** variance: the Poisson signal part of each $\sigma^{\mathrm{corr}}_{ij}$ is rebuilt at the reflection's merged $\langle I\rangle$ rather than at that observation's own intensity. Weighting by an observation's own $\sigma^2$ biases the inverse-variance mean low below about one photon, because an up-fluctuated observation gets a larger sigma and is then down-weighted too hard. The rotation combine already does this; for stills it is on by default, and `--no-expected-variance-merge` restores the observed-sigma weighting. + An internal-consistency term can inflate uncertainties when multiple observations are present, in the spirit of XSCALE. ### 10.5 Merging statistics @@ -600,7 +604,7 @@ After scale-fulls, three **correction surfaces** are fitted on the combined full - **Decay.** Radiation damage weakens later frames more at higher resolution — a resolution×time (Debye–Waller) systematic the resolution-flat per-image scale cannot capture. A single global relative-$B$ rate is fitted, $\ln(I_\mathrm{ref}/I_\mathrm{obs}) = 2\,(\mathrm{d}B/\mathrm{d}n)\,(n-\bar n)\,s^2$ (frame $n$, $s^2 = 1/4d^2$), and folded into the scale. It engages only when the total relative-$B$ over the run exceeds a physical floor (2 Ų); below that the decay is negligible and "correcting" it only spreads symmetry equivalents (same $s^2$, different frames). An optional **per-batch relative-$B$** (`--relative-b[=deg]`, off unless requested; 10°-of-rotation batches by default) extends the single global rate to a smooth $B(n)$ curve — the same $s^2$-weighted decay fit solved independently over short frame batches, curvature-penalized so it cannot over-fit and cross-validated like the surfaces below — for crystals whose decay is non-linear in dose. - **Absorption.** A smooth multiplicative factor over the diffracted-beam direction expressed in the goniometer (crystal) frame: each full's predicted detector position gives the lab diffracted direction, de-rotated by the spindle so a fixed crystal-frame direction is sampled at many rotation angles and its grid cell is well-determined. Negligible at hard X-rays / thin crystals; it matters at low photon energy. -- **Modulation** (detector-plane flat-field). A smooth multiplicative factor over where each reflection lands on the detector (predicted $x,y$): symmetry-equivalents land at different positions as the crystal rotates, over-determining the surface. It absorbs detector-response and geometric systematics that inflate $R_\mathrm{meas}$. The same 16×16 detector-frame surface is available for the stills path (`--stills-modulation`, off by default), where serial data repeatedly hammers the same detector regions. +- **Modulation** (detector-plane flat-field). A smooth multiplicative factor over where each reflection lands on the detector (predicted $x,y$): symmetry-equivalents land at different positions as the crystal rotates, over-determining the surface. It absorbs detector-response and geometric systematics that inflate $R_\mathrm{meas}$. Each surface is **cross-validated**: fitted on even-numbered frames and kept only if it improves the held-out odd-frame agreement by a clear margin (and vice versa), scored by a **σ-independent, $R_\mathrm{meas}$-like** fractional agreement $\sum|I_s-I_\mathrm{ref}|/\sum|I_\mathrm{ref}|$ rather than a studentized $\chi^2$ — so a surface cannot "pass" by reshaping the sigmas instead of tightening the intensities. A surface fitted to noise where its systematic is absent does not generalize and is discarded — a correction never adds scatter. @@ -695,10 +699,10 @@ A **dataset-wide** Wilson $B$ is also estimated over the merged reflections — - **Bragg integration is profile-fitted by default** (per-shell Gaussian profile, Kabsch extraction; §9.3), with plain box summation available as a fallback (`--integrator boxsum`). The profiles are built per frame from that frame's strong spots, which suits fast-feedback and serial/streaming use; a profile shared across many frames (as in full offline workflows) is not currently formed. - **Space-group symmetry** beyond centering absences is not necessarily enforced during prediction/integration unless the space group is supplied and used downstream. - **Resolution masking and ice rings** are controllable; including ice-ring spots in indexing can improve robustness for some samples but may bias refinement in others. -- **Rotation vs still modes** differ substantially in prediction and scaling: partiality is angle-driven in rotation data, while stills are predicted (within an excitation-error window) and scaled with unit partiality. +- **Rotation vs still modes** differ substantially in prediction and scaling: partiality is angle-driven in rotation data, while stills are predicted within an excitation-error window and get their partiality from the default-on per-crystal tilt post-refinement (§10.2) — or unit partiality with `--simple-stills`. - **Space-group determination.** When no space group is supplied, a POINTLESS-like search scores Laue-group symmetry (CC of $I(h)$ vs $I(Rh)$ plus merge self-consistency) and detects screw/centering absences from the $P1$-merged intensities. The self-consistency test is calibrated so a merohedral twin — whose twin law forces non-equivalent reflections together and inflates the merged $\chi^2$ — stays in its true lower symmetry rather than being over-promoted to the holohedral group. Because a partial twin's within-orbit $\chi^2$ can nonetheless look self-consistent, a chi²-passing promotion is additionally **vetoed** when merging its extra operator balloons the error-model $b$ (the intensity-proportional systematic) relative to the confirmed subgroup: a genuine symmetry step gains multiplicity without inflating $b$, whereas a twin forces non-equivalent reflections together and $b$ balloons. **Centering** is accepted when the systematically-absent class is weak relative to the present one by *either* of two floor-independent tests — its mean signed $I/\sigma$ well below the present mean, *or* its rate of individually-significant reflections well below the present class's own significant rate. The second test matters on weak / low-energy data, where a positive intensity floor (background/profile leakage) lifts the absent class's mean $I/\sigma$ to $\sim1.5$–$2.3$ instead of $\sim0$ and, when the present class is itself weak, inflates the plain mean ratio past its bound and hides a real centering (an $I$-centred cubic recorded at 5 keV was otherwise kept primitive); a false centering fails both tests because its absent class is as strong as the present one. When several centerings pass, they are ranked by their **net** systematic absences (absent minus violating), not the gross absent count, so a super-centering (e.g. $F$ over a true $C$) whose extra, only-half-populated absent class merely dilutes the strength ratio does not out-rank the correct lower centering. - **Twinning check.** A Padilla–Yeates $L$-test ($\langle|L|\rangle$, $\langle L^2\rangle$) and the second moment $\langle I^2\rangle/\langle I\rangle^2$ (taken per resolution shell with noise-only shells skipped and Wilson outliers rejected, so a single strong reflection in a collapsed-mean shell cannot skew it) are written to the merged mmCIF as a twinning diagnostic. Twinning is only flagged in Laue classes where a merohedral twin law can exist; the holohedral high-symmetry classes ($4/mmm$, $6/mmm$, $m\bar{3}m$, and $\bar{3}m$ on a rhombohedral lattice) are exempt, so a low $\langle|L|\rangle$ there is reported as a statistical artefact rather than twinning. -- **Outlier rejection.** Merging applies an optional per-observation median-based $N\sigma$ cut (default 6σ for `rot3d`) and an optional per-crystal $\Delta\mathrm{CC}_{1/2}$ image rejection (`--reject-delta-cchalf`, CrystFEL-style, off by default). The same $N\sigma$ cut is fed back into the error model: after an initial $a,b$ fit the parameters are re-fit once on the reflections that survive rejection (dropping any whose squared deviation exceeds $N\sigma^2\,[a\,\sigma^2 + (b\,\langle I\rangle)^2]$), so the calibrated errors describe the reflections that actually enter the merge rather than the pre-rejection pool. +- **Outlier rejection.** Merging applies an optional per-observation median-based $N\sigma$ cut (`--reject-outliers`, default 6σ for `rot3d`, off otherwise). The same $N\sigma$ cut is fed back into the error model: after an initial $a,b$ fit the parameters are re-fit once on the reflections that survive rejection (dropping any whose squared deviation exceeds $N\sigma^2\,[a\,\sigma^2 + (b\,\langle I\rangle)^2]$), so the calibrated errors describe the reflections that actually enter the merge rather than the pre-rejection pool. - **Automatic resolution cutoff.** By default the reported/written high-resolution limit is trimmed where $\mathrm{CC}_{1/2}$ falls off (logistic, target 0.30); `--scaling-high-resolution` overrides it and `--resolution-cutoff off` disables it. - **Amplitudes and intensities.** The merged output carries both intensities (mmCIF `intensity_meas`, MTZ `IMEAN`/`SIGIMEAN`) and French–Wilson amplitudes (mmCIF `F_meas_au`, MTZ `F`/`SIGF`; §10.8), so a downstream program can refine against either. diff --git a/docs/HDF5.md b/docs/HDF5.md index e561730e..2d095956 100644 --- a/docs/HDF5.md +++ b/docs/HDF5.md @@ -261,7 +261,6 @@ In legacy/VDS mode these live in the data files and are linked/virtual-stacked i | `imageScaleFactor` | | on-the-fly per-image scale factor *g* | | `imageScaleCC` | | on-the-fly scaling correlation coefficient | | `imageScaleMosaicity` | deg | scaling-model mosaicity | -| `imageScaleBFactor` | Ų | scaling-model B-factor | **Per-image lattices:** `latticeIndexed` `[n_images, 9]` (Å) — the real-space lattice (flattened 3×3); `latticeIndexedExtra` `[n_images, max_extra_lattices, 9]` (Å) — additional orientation diff --git a/docs/RUGNUX.md b/docs/RUGNUX.md index 68075a24..5bb47bff 100644 --- a/docs/RUGNUX.md +++ b/docs/RUGNUX.md @@ -150,14 +150,13 @@ 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 \ - --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 -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 +`ffbidx` requires a known cell (`-C`) and is the indexer of choice for sparse serial stills. The +self-calibrating spot finder is on by default for stills (`--no-adaptive-spots` turns it off), and for +serial stills 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 @@ -189,9 +188,10 @@ 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) | +| `--adaptive-spots` | Self-calibrating detection (**default for stills**): the strong-pixel threshold comes from each image's own per-resolution-ring noise instead of the fixed `--spot-threshold`, so one setting adapts across datasets (no per-dataset `--spot-threshold`/`--spot-sigma` tuning). Rotation data keeps the fixed-threshold finder unless this is given | +| `--no-adaptive-spots` | Turn adaptive detection off and use the fixed `--spot-threshold` / `--spot-sigma` finder | | `--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-high-resolution ` | High-resolution limit for spot finding, Å. Omitted: stills extend as far as the detector reaches (no resolution clipping); rotation data keeps a 1.5 Å limit | | `--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. **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) | @@ -203,7 +203,7 @@ Azimuthal integration (the radial profile behind the per-image ice-ring score): | --- | --- | | `-q, --azim-q-spacing ` | Q bin spacing, 1/Å (default: 0.01; finer resolves the narrow ice rings) | | `--azim-min-q ` | Minimum Q, 1/Å | -| `--azim-max-q ` | Maximum Q, 1/Å | +| `--azim-max-q ` | Maximum Q, 1/Å. Omitted: integration extends to the highest Q the detector reaches. The adaptive spot finder shares these Q bins, so this also sets how far self-calibrating detection can see | | `--azim-phi-bins ` | Number of azimuthal (phi) bins (default: 1) | | `--polarization-correction ` | Enable/disable the azimuthal polarization correction | | `--solid-angle-correction ` | Enable/disable the azimuthal solid-angle correction | @@ -239,14 +239,12 @@ Scaling and merging: | --- | --- | | `--no-merge` | Skip scaling and merging (on by default); write only the per-image `_process.h5` | | `-A, --anomalous` | Anomalous mode (keep Friedel pairs separate) | -| `-B, --refine-bfactor` | Refine a per-image B-factor (stills only) | | `--scale-fulls` / `--no-scale-fulls` | rot3d: refit a per-frame scale on the combined fulls (XDS order, Unity model); on by default for rotation data, off for stills | | `--smooth-g[=deg]` | rot3d: smooth the per-frame scale *G* over a degree range before the 3D combine (XDS DELPHI-like; default 5° for rotation, 0 = off) | | `--no-scaling-corrections` | rot3d: disable the default-on decay + absorption + modulation correction surfaces fitted on the fulls after scale-fulls (see below) | | `--relative-b[=deg]` | rot3d: fit a per-batch relative-*B* beyond the single decay slope over deg-degree batches, cross-validated (default 10° when bare; off otherwise) | -| `--still-partiality` | Experimental (stills): weight reflections by a Gaussian excitation-error partiality instead of treating each as a full | -| `--partiality-uncertainty ` | Stills: extra merge sigma ~num·(1−partiality)·⟨I⟩ on partials (use with `--still-partiality`; default 0, ~2.5 recommended) | -| `--stills-modulation` | Experimental (stills): fit a detector-plane modulation (flat-field) surface, cross-validated (default off) | +| `--simple-stills` | Stills: treat every reflection as a full (*p* = 1, single-pass scale/merge) — disables the default-on physical partiality post-refinement | +| `--no-expected-variance-merge` | Stills: disable the default expected-variance merge weighting (which rebuilds each weak observation's signal variance at the reflection mean to de-bias the inverse-variance merge); restores observed-sigma weighting | | `--capture-uncertainty ` | rot3d: systematic sigma on under-captured fulls, ~num·(1−captured_fraction)·I (default: 1.0 for rotation, 0 otherwise) | | `--min-captured-fraction ` | rot3d: drop a combined full whose rocking curve was captured below this fraction — edge-of-sweep truncated fulls (default: 0.7 for rotation, 0 otherwise; 0 = off) | | `--scaling-high-resolution ` | High-resolution limit for scaling, Å — manual override (default: no limit; disables the automatic cutoff below) | @@ -255,7 +253,6 @@ Scaling and merging: | `--resolution-shells ` | Number of resolution shells in the reported statistics table (default: 10) | | `--min-partiality ` | Minimum partiality to accept a reflection (default: 0.02) | | `--reject-outliers ` | Per-observation outlier rejection, N σ from the per-reflection median (default: 6 for `rot3d`, off otherwise) | -| `--reject-delta-cchalf ` | Drop images with ΔCC1/2 below mean − N·stddev (default: off) | | `--min-image-cc ` | Per-image CC limit, percent (default: no limit) | | `--mosaicity ` | Diagnostic: fix the scaling mosaicity (°) instead of using the per-image seed | | `--scaling-iterations ` | Scaling iterations with no reference data (default: 3) | diff --git a/fpga/pcie_driver/postinstall.sh b/fpga/pcie_driver/postinstall.sh index 6340c150..2ce71017 100644 --- a/fpga/pcie_driver/postinstall.sh +++ b/fpga/pcie_driver/postinstall.sh @@ -5,7 +5,7 @@ VERSION=1.0.0-rc.160 occurrences=`/usr/sbin/dkms status | grep jfjoch | grep ${VERSION} | wc -l` -if [ ! occurrences > 0 ]; then +if [ "$occurrences" -eq 0 ]; then /usr/sbin/dkms add -m jfjoch -v ${VERSION} fi /usr/sbin/dkms build -m jfjoch -v ${VERSION} diff --git a/frame_serialize/CBORStream2Deserializer.cpp b/frame_serialize/CBORStream2Deserializer.cpp index 9c36c8a8..14642a35 100644 --- a/frame_serialize/CBORStream2Deserializer.cpp +++ b/frame_serialize/CBORStream2Deserializer.cpp @@ -828,8 +828,6 @@ namespace { message.image_scale_cc = GetCBORFloat(value); else if (key == "image_scale_mosaicity") message.image_scale_mosaicity = GetCBORFloat(value); - else if (key == "image_scale_b_factor") - message.image_scale_b_factor = GetCBORFloat(value); else if (key == "roi_integrals") ProcessROIElementMap(message, value); else { @@ -1449,8 +1447,6 @@ namespace { GetCBORFloatArray(value, message.image_scale_cc); else if (key == "image_scale_mosaicity") GetCBORFloatArray(value, message.image_scale_mosaicity); - else if (key == "image_scale_b_factor") - GetCBORFloatArray(value, message.image_scale_b_factor); else if (key == "integrated_reflections") GetCBORInt32Array(value, message.integrated_reflections); else if (key == "niggli_class") diff --git a/frame_serialize/CBORStream2Serializer.cpp b/frame_serialize/CBORStream2Serializer.cpp index da5134fe..813bedf7 100644 --- a/frame_serialize/CBORStream2Serializer.cpp +++ b/frame_serialize/CBORStream2Serializer.cpp @@ -785,7 +785,6 @@ void CBORStream2Serializer::SerializeSequenceEnd(const EndMessage& message) { CBOR_ENC(mapEncoder, "image_scale_factor", message.image_scale_factor); CBOR_ENC(mapEncoder, "image_scale_cc", message.image_scale_cc); CBOR_ENC(mapEncoder, "image_scale_mosaicity", message.image_scale_mosaicity); - CBOR_ENC(mapEncoder, "image_scale_b_factor", message.image_scale_b_factor); CBOR_ENC(mapEncoder, "integrated_reflections", message.integrated_reflections); CBOR_ENC(mapEncoder, "niggli_class", message.niggli_class); CBOR_ENC(mapEncoder, "pixel_sum", message.pixel_sum); @@ -870,7 +869,6 @@ void CBORStream2Serializer::SerializeImageInternal(CborEncoder &mapEncoder, cons CBOR_ENC(mapEncoder, "beam_corr_y", message.beam_corr_y); CBOR_ENC(mapEncoder, "image_scale_factor", message.image_scale_factor); CBOR_ENC(mapEncoder, "image_scale_mosaicity", message.image_scale_mosaicity); - CBOR_ENC(mapEncoder, "image_scale_b_factor", message.image_scale_b_factor); CBOR_ENC(mapEncoder, "image_scale_cc", message.image_scale_cc); CBOR_ENC(mapEncoder, "user_data", message.user_data.dump()); diff --git a/frontend/package-lock.json b/frontend/package-lock.json index 02bf59d5..08759dc6 100644 --- a/frontend/package-lock.json +++ b/frontend/package-lock.json @@ -1,12 +1,12 @@ { "name": "jungfraujoch-frontend", - "version": "1.0.0-rc.159", + "version": "1.0.0-rc.160", "lockfileVersion": 3, "requires": true, "packages": { "": { "name": "jungfraujoch-frontend", - "version": "1.0.0-rc.159", + "version": "1.0.0-rc.160", "license": "GPL-3.0", "dependencies": { "@emotion/react": "^11.10.4", diff --git a/frontend/src/client/types.gen.ts b/frontend/src/client/types.gen.ts index eba0e361..c2095aa0 100644 --- a/frontend/src/client/types.gen.ts +++ b/frontend/src/client/types.gen.ts @@ -712,9 +712,11 @@ export type spot_finding_settings = { min_pix_per_spot: number; max_pix_per_spot: number; /** - * High resolution limit for spot finding [Angstrom] + * High resolution limit for spot finding [Angstrom]. Optional: if omitted, spot finding extends + * as far as the detector reaches, i.e. the detection is not clipped in resolution. + * */ - high_resolution_limit: number; + high_resolution_limit?: number; /** * Low resolution limit for spot finding [Angstrom] */ @@ -752,7 +754,13 @@ export type azim_int_settings = { * Apply solid angle correction for azimuthal integration */ solid_angle_corr: boolean; - high_q_recipA: number; + /** + * Upper q limit of the azimuthal integration [1/Angstrom]. Optional: if omitted, the integration + * (and the adaptive spot detection that shares these q bins) extends to the highest q the + * detector reaches. + * + */ + high_q_recipA?: number; low_q_recipA: number; q_spacing: number; /** @@ -1741,7 +1749,6 @@ export const plot_type = { INTEGRATED_REFLECTIONS: 'integrated_reflections', IMAGE_SCALE_FACTOR: 'image_scale_factor', IMAGE_SCALE_CC: 'image_scale_cc', - IMAGE_SCALE_B: 'image_scale_b', COMPRESSION_RATIO: 'compression_ratio', ICE_RING_SCORE: 'ice_ring_score' } as const; @@ -3043,7 +3050,7 @@ export type getPreviewPlotData = { /** * Type of requested plot */ - type: 'bkg_estimate' | 'azint' | 'azint_1d' | 'spot_count' | 'spot_count_low_res' | 'spot_count_indexed' | 'spot_count_ice' | 'indexing_rate' | 'indexing_lattice_count' | 'indexing_unit_cell_length' | 'indexing_unit_cell_angle' | 'profile_radius' | 'mosaicity' | 'b_factor' | 'error_pixels' | 'saturated_pixels' | 'image_collection_efficiency' | 'receiver_delay' | 'receiver_free_send_buf' | 'strong_pixels' | 'roi_sum' | 'roi_mean' | 'roi_max_count' | 'roi_pixels' | 'roi_weighted_x' | 'roi_weighted_y' | 'packets_received' | 'max_pixel_value' | 'resolution_estimate' | 'pixel_sum' | 'processing_time' | 'beam_center_x' | 'beam_center_y' | 'integrated_reflections' | 'image_scale_factor' | 'image_scale_cc' | 'image_scale_b' | 'compression_ratio' | 'ice_ring_score'; + type: 'bkg_estimate' | 'azint' | 'azint_1d' | 'spot_count' | 'spot_count_low_res' | 'spot_count_indexed' | 'spot_count_ice' | 'indexing_rate' | 'indexing_lattice_count' | 'indexing_unit_cell_length' | 'indexing_unit_cell_angle' | 'profile_radius' | 'mosaicity' | 'b_factor' | 'error_pixels' | 'saturated_pixels' | 'image_collection_efficiency' | 'receiver_delay' | 'receiver_free_send_buf' | 'strong_pixels' | 'roi_sum' | 'roi_mean' | 'roi_max_count' | 'roi_pixels' | 'roi_weighted_x' | 'roi_weighted_y' | 'packets_received' | 'max_pixel_value' | 'resolution_estimate' | 'pixel_sum' | 'processing_time' | 'beam_center_x' | 'beam_center_y' | 'integrated_reflections' | 'image_scale_factor' | 'image_scale_cc' | 'compression_ratio' | 'ice_ring_score'; /** * Fill value for elements that were missed during data collection * @@ -3090,7 +3097,7 @@ export type getPreviewPlotBinData = { /** * Type of requested plot */ - type: 'bkg_estimate' | 'azint' | 'azint_1d' | 'spot_count' | 'spot_count_low_res' | 'spot_count_indexed' | 'spot_count_ice' | 'indexing_rate' | 'indexing_lattice_count' | 'indexing_unit_cell_length' | 'indexing_unit_cell_angle' | 'profile_radius' | 'mosaicity' | 'b_factor' | 'error_pixels' | 'saturated_pixels' | 'image_collection_efficiency' | 'receiver_delay' | 'receiver_free_send_buf' | 'strong_pixels' | 'roi_sum' | 'roi_mean' | 'roi_max_count' | 'roi_pixels' | 'roi_weighted_x' | 'roi_weighted_y' | 'packets_received' | 'max_pixel_value' | 'resolution_estimate' | 'pixel_sum' | 'processing_time' | 'beam_center_x' | 'beam_center_y' | 'integrated_reflections' | 'image_scale_factor' | 'image_scale_cc' | 'image_scale_b' | 'compression_ratio' | 'ice_ring_score'; + type: 'bkg_estimate' | 'azint' | 'azint_1d' | 'spot_count' | 'spot_count_low_res' | 'spot_count_indexed' | 'spot_count_ice' | 'indexing_rate' | 'indexing_lattice_count' | 'indexing_unit_cell_length' | 'indexing_unit_cell_angle' | 'profile_radius' | 'mosaicity' | 'b_factor' | 'error_pixels' | 'saturated_pixels' | 'image_collection_efficiency' | 'receiver_delay' | 'receiver_free_send_buf' | 'strong_pixels' | 'roi_sum' | 'roi_mean' | 'roi_max_count' | 'roi_pixels' | 'roi_weighted_x' | 'roi_weighted_y' | 'packets_received' | 'max_pixel_value' | 'resolution_estimate' | 'pixel_sum' | 'processing_time' | 'beam_center_x' | 'beam_center_y' | 'integrated_reflections' | 'image_scale_factor' | 'image_scale_cc' | 'compression_ratio' | 'ice_ring_score'; /** * Name of ROI for which plot is requested */ diff --git a/frontend/src/client/zod.gen.ts b/frontend/src/client/zod.gen.ts index e74daad1..2a71e523 100644 --- a/frontend/src/client/zod.gen.ts +++ b/frontend/src/client/zod.gen.ts @@ -300,7 +300,7 @@ export const zSpotFindingSettings = z.object({ photon_count_threshold: z.coerce.bigint().gte(BigInt(0)).max(BigInt('9223372036854775807'), { error: 'Invalid value: Expected int64 to be <= 9223372036854775807' }), min_pix_per_spot: z.coerce.bigint().gte(BigInt(1)).max(BigInt('9223372036854775807'), { error: 'Invalid value: Expected int64 to be <= 9223372036854775807' }), max_pix_per_spot: z.coerce.bigint().gte(BigInt(1)).max(BigInt('9223372036854775807'), { error: 'Invalid value: Expected int64 to be <= 9223372036854775807' }), - high_resolution_limit: z.number(), + high_resolution_limit: z.number().optional(), low_resolution_limit: z.number(), high_resolution_limit_for_spot_count_low_res: z.number().gte(2).lte(8), quick_integration: z.boolean().default(false), @@ -311,7 +311,7 @@ export const zSpotFindingSettings = z.object({ export const zAzimIntSettings = z.object({ polarization_corr: z.boolean().default(true), solid_angle_corr: z.boolean().default(true), - high_q_recipA: z.number().gte(0.00002).lte(10), + high_q_recipA: z.number().gte(0.00002).lte(10).optional(), low_q_recipA: z.number().gte(0.00001).lte(10), q_spacing: z.number().gte(0.00001), azimuthal_bins: z.coerce.bigint().gte(BigInt(1)).lte(BigInt(512)).optional().default(BigInt(1)), @@ -840,7 +840,6 @@ export const zPlotType = z.enum([ 'integrated_reflections', 'image_scale_factor', 'image_scale_cc', - 'image_scale_b', 'compression_ratio', 'ice_ring_score' ]); @@ -1150,7 +1149,6 @@ export const zGetPreviewPlotQuery = z.object({ 'integrated_reflections', 'image_scale_factor', 'image_scale_cc', - 'image_scale_b', 'compression_ratio', 'ice_ring_score' ]), @@ -1206,7 +1204,6 @@ export const zGetPreviewPlotBinQuery = z.object({ 'integrated_reflections', 'image_scale_factor', 'image_scale_cc', - 'image_scale_b', 'compression_ratio', 'ice_ring_score' ]), diff --git a/frontend/src/components/AzIntSettings.tsx b/frontend/src/components/AzIntSettings.tsx index c84ca9fb..73beedd9 100644 --- a/frontend/src/components/AzIntSettings.tsx +++ b/frontend/src/components/AzIntSettings.tsx @@ -30,6 +30,7 @@ function AzIntSettings({s: serverS}: MyProps) { const [lowQError, setLowQError] = useState(false); const [highQError, setHighQError] = useState(false); const [qSpacingError, setQSpacingError] = useState(false); + const [highQ, setHighQ] = useState(5); const { submit, pending, snackbar } = useUpload(putConfigAzimIntMutation()); useEffect(() => { @@ -37,6 +38,7 @@ function AzIntSettings({s: serverS}: MyProps) { setS(serverS); setLastDownloadedS(serverS); setDownloadCounter(c => c + 1); + setHighQ(serverS.high_q_recipA ?? highQ); } // eslint-disable-next-line react-hooks/exhaustive-deps }, [serverS]); @@ -46,7 +48,8 @@ function AzIntSettings({s: serverS}: MyProps) { return ( submit(s)} - uploadDisabled={pending || (s.high_q_recipA <= s.low_q_recipA) || highQError || lowQError || qSpacingError} + uploadDisabled={pending || ((s.high_q_recipA !== undefined) && (s.high_q_recipA <= s.low_q_recipA)) + || highQError || lowQError || qSpacingError} snackbar={snackbar}> { + setHighQ(val); setS(prev => ({...prev, high_q_recipA: val})); setHighQError(err); }} fullWidth/> + ) => { + // No high Q means integrating out to the highest Q the detector reaches. + setS(prev => event.target.checked ? _.omit(prev, "high_q_recipA") + : {...prev, high_q_recipA: highQ}); + }} + />} label="High Q to detector edge"/> (default_spot_finding_settings); const [lastDownloadedS, setLastDownloadedS] = useState(default_spot_finding_settings); const [highResGap, setHighResGap] = useState(1.5); + const [highResLimit, setHighResLimit] = useState(2.5); // Only adopt the server copy when it actually changed, otherwise the // 1 s statistics poll would overwrite edits the user is making. @@ -42,6 +43,7 @@ function DataProcessingSettings({s: serverS, update}: MyProps) { setS(serverS); setLastDownloadedS(serverS); setHighResGap(serverS.high_res_gap_Q_recipA ?? highResGap); + setHighResLimit(serverS.high_resolution_limit ?? highResLimit); } // eslint-disable-next-line react-hooks/exhaustive-deps }, [serverS]); @@ -64,8 +66,15 @@ function DataProcessingSettings({s: serverS, update}: MyProps) { const setLowResolutionLimit = (event: Event, newValue: number | number[]) => apply({...s, low_resolution_limit: newValue as number}); - const setHighResolutionLimit = (event: Event, newValue: number | number[]) => - apply({...s, high_resolution_limit: newValue as number}); + const setHighResolutionLimit = (event: Event, newValue: number | number[]) => { + const v = newValue as number; + setHighResLimit(v); + apply({...s, high_resolution_limit: v}); + }; + + // No high-resolution limit means spot finding goes as far as the detector reaches. + const autoHighResolutionLimitToggle = (event: ChangeEvent) => + apply({...s, high_resolution_limit: event.target.checked ? undefined : highResLimit}); const setHighResolutionLimitForCountingLowResSpots = (event: Event, newValue: number | number[]) => apply({...s, high_resolution_limit_for_spot_count_low_res: newValue as number}); @@ -131,9 +140,13 @@ function DataProcessingSettings({s: serverS, update}: MyProps) { valueLabelFormat={(value) => value.toFixed(1)} /> + + To detector edge High resolution limit [Å] - value.toFixed(1)} diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 08956bec..3e5e6885 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -104,6 +104,7 @@ IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetry(Data // already ordered non-ice-first, strongest-first (FilterSpotsByCount), so the prefix IS the seed. const float idx_tol = experiment.GetIndexingSettings().GetTolerance(); const float idx_tol_sq = idx_tol * idx_tol; + constexpr float SEED_STOP_FRACTION = 0.9f; // seed explained this well -> stop escalating IndexerResult indexer_result; bool any_executed = false; float best_frac = -1.0f; @@ -133,6 +134,11 @@ IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetry(Data } const float frac = recip.empty() ? 0.0f : static_cast(n) / recip.size(); if (frac > best_frac) { best_frac = frac; indexer_result = std::move(res); } + // A lattice that already explains nearly the whole seed is kept whatever a larger seed + // returns: the winner is the highest explained FRACTION, and adding weaker spots almost + // always lowers it. Stop here - this is what keeps clean frames at one indexer call. + if (frac >= SEED_STOP_FRACTION) + break; } if (recip.size() < seed_cap) // already fed every available spot; a larger cap won't add any break; @@ -408,7 +414,6 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg, .geom = outcome.experiment.GetDiffractionGeometry(), .latt = latt, .mosaicity_deg = mos_deg, - .image_scale_b_factor_Ang2 = msg.image_scale_b_factor, .image_scale_cc = msg.image_scale_cc, }; @@ -425,10 +430,6 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg, .mosaicity_deg = std::fabs(mos_deg), // FWHM -> sigma; 0 when monochromatic, leaving the prediction unchanged. .bandwidth_sigma = experiment.GetBandwidthFWHM().value_or(0.0f) / 2.3548f, - // Experimental stills partiality (off by default): sigma = ewald_dist_cutoff/2 = the per-image - // profile radius, so a reflection at the acceptance edge (dist_ewald ~ 2*sigma) keeps p ~ exp(-2). - .still_partiality = experiment.GetBraggIntegrationSettings().GetStillPartiality(), - .profile_radius_recipA = ewald_dist_cutoff * 0.5f }; // Predict, then integrate with the selected integrator (box-sum or profile-fit). @@ -510,14 +511,13 @@ IndexAndRefine::DetermineRefineAnalyze(DataMessage &msg, const SpotFindingSettin return outcome; } -bool IndexAndRefine::ProcessImage(DataMessage &msg, +void 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 dda9a0d3..54d09b51 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -104,7 +104,7 @@ public: // 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, + void 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/IntegrationOutcome.h b/image_analysis/IntegrationOutcome.h index 3b9cbec6..75b2a09c 100644 --- a/image_analysis/IntegrationOutcome.h +++ b/image_analysis/IntegrationOutcome.h @@ -12,7 +12,6 @@ struct IntegrationOutcome { CrystalLattice latt; std::vector reflections; std::optional mosaicity_deg; - std::optional image_scale_b_factor_Ang2; std::optional image_scale_cc; std::optional image_scale_cc_n; std::optional image_scale_g; diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index 7d188a38..5070dbf8 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -134,35 +134,53 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, && 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(); + // Choose the per-image min-pix adaptively instead of a fixed one. min-pix filters connected + // components AFTER detection, so detection (the expensive per-pixel pass) runs ONCE and only + // the cheap CCL + spot filter is repeated; the azimuthal profile is the one Detect() computed. + // Index at 3/2/1 (index-only, no integration/accumulation), keep whichever maximises + // n_indexed^2 / n_total (indexed count weighted by indexed fraction) together with its spot + // list, and integrate that one. Keeping the list also means the frame integrated is exactly + // the frame scored, which re-extracting could not guarantee on the GPU (atomic-order sums). + const auto detect_start_time = std::chrono::steady_clock::now(); + finder.Detect(*preprocessor_buffer, spot_finding_settings); + float spot_finding_time_s = + std::chrono::duration(std::chrono::steady_clock::now() - detect_start_time).count(); + float indexing_time_s = 0.0f; + SpotFindingSettings s = spot_finding_settings; + std::vector best_spots; 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); + const auto extract_start_time = std::chrono::steady_clock::now(); + std::vector spots = finder.ExtractSpots(*preprocessor_buffer, s, mask_resolution); + spot_finding_time_s += + std::chrono::duration(std::chrono::steady_clock::now() - extract_start_time).count(); SpotAnalyze(experiment, s, spots, output); - if (indexer.IndexFrameOnly(output, s)) { + const bool indexed = indexer.IndexFrameOnly(output, s); + indexing_time_s += output.indexing_time_s.value_or(0.0f); + if (indexed) { 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 (score > best_score) { + best_score = score; + best_mp = mp; + best_spots = std::move(spots); + } } } if (best_mp != 0) { - // Re-run spot finding + index at the winning min-pix and integrate there. + // Index and integrate the winning spot list; no spot finding left to do. s.min_pix_per_spot = best_mp; - const std::vector spots = finder.Run(*preprocessor_buffer, s, mask_resolution); - SpotAnalyze(experiment, s, spots, output); + SpotAnalyze(experiment, s, best_spots, output); indexer.ProcessImage(output, s, *prediction, integrate_fn); + indexing_time_s += output.indexing_time_s.value_or(0.0f); } - output.spot_finding_time_s = std::chrono::duration(std::chrono::steady_clock::now() - start_time).count(); + // Each indexer call reports only its own time, so the escalation's total is summed here. + output.spot_finding_time_s = spot_finding_time_s; + output.indexing_time_s = indexing_time_s; } 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); @@ -174,8 +192,8 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, #ifdef JFJOCH_USE_CUDA if (fused) { - // 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. + // Lift the azimuthal profile the fused engine computed in the same detection pass; its azint + // cost is folded into spot_finding_time_s above. profile.Clear(integration); profile += fused_adaptive->GetProfile(); output.azint_time_s = 0.0f; @@ -230,7 +248,10 @@ void MXAnalysisWithoutFPGA::RunROIOnly(DataMessage &output) { void MXAnalysisWithoutFPGA::UpdateMaskResolution(const SpotFindingSettings &settings) { mask_low_res = settings.low_resolution_limit; mask_high_res = settings.high_resolution_limit; + // No high-resolution limit requested -> mask nothing at the high-resolution end: no pixel has d < 0, + // and the detector's own edge is where the pixels stop anyway. + const float high_res = mask_high_res.value_or(0.0f); auto const &resolution_map = integration.Resolution(); for (int i = 0; i < mask_resolution.size(); i++) - mask_resolution[i] = (resolution_map[i] > mask_low_res) || (resolution_map[i] < mask_high_res); + mask_resolution[i] = (resolution_map[i] > mask_low_res) || (resolution_map[i] < high_res); } diff --git a/image_analysis/MXAnalysisWithoutFPGA.h b/image_analysis/MXAnalysisWithoutFPGA.h index c79f473a..64eb08f7 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.h +++ b/image_analysis/MXAnalysisWithoutFPGA.h @@ -54,7 +54,9 @@ class MXAnalysisWithoutFPGA { const PixelMask &mask; std::vector mask_resolution; - float mask_high_res; + // The limits mask_resolution was built for. Kept as the OPTIONAL the caller passed, so an unset + // high-resolution limit compares equal to itself and the mask is not rebuilt on every image. + std::optional mask_high_res; float mask_low_res; void UpdateMaskResolution(const SpotFindingSettings& settings); #ifdef JFJOCH_USE_CUDA diff --git a/image_analysis/bragg_prediction/BraggPrediction.cpp b/image_analysis/bragg_prediction/BraggPrediction.cpp index e19d0249..ac483c60 100644 --- a/image_analysis/bragg_prediction/BraggPrediction.cpp +++ b/image_analysis/bragg_prediction/BraggPrediction.cpp @@ -156,13 +156,6 @@ int BraggPrediction::Calc(const DiffractionExperiment &experiment, const Crystal continue; float d = 1.0f / sqrtf(recip_sq); - float partiality = 1.0f; - if (settings.still_partiality && settings.profile_radius_recipA > 0.0f) { - const float sig_bw = settings.bandwidth_sigma * std::fabs(recip_z); - const float sigma2 = settings.profile_radius_recipA * settings.profile_radius_recipA - + sig_bw * sig_bw; - partiality = std::exp(-0.5f * dist_ewald_sphere * dist_ewald_sphere / sigma2); - } reflections[i] = Reflection{ .h = h, .k = k, @@ -175,7 +168,7 @@ int BraggPrediction::Calc(const DiffractionExperiment &experiment, const Crystal .d = d, .dist_ewald = dist_ewald_sphere, .rlp = 1.0, - .partiality = partiality, + .partiality = 1.0f, .zeta = 1.0, .image_scale_corr = 1.0 }; diff --git a/image_analysis/bragg_prediction/BraggPrediction.h b/image_analysis/bragg_prediction/BraggPrediction.h index 5406d23e..60400b00 100644 --- a/image_analysis/bragg_prediction/BraggPrediction.h +++ b/image_analysis/bragg_prediction/BraggPrediction.h @@ -23,13 +23,6 @@ struct BraggPredictionSettings { // σ_bw = |recip_z|·bandwidth_sigma (= bλ/2d²), so the 1/d² pink-beam smear no // longer clips high-resolution reflections. float bandwidth_sigma = 0.0f; - // Experimental stills partiality (rugnux --still-partiality). When still_partiality is set and - // profile_radius_recipA > 0, each reflection gets a Gaussian excitation-error partiality - // p = exp(-dist_ewald^2 / (2*sigma^2)), sigma^2 = profile_radius_recipA^2 + (bandwidth_sigma*|recip_z|)^2, - // instead of the fixed 1.0 (full). Off by default. Keep these two as the trailing members so the - // existing designated initializers (which stop at bandwidth_sigma) remain valid. - bool still_partiality = false; - float profile_radius_recipA = 0.0f; }; class BraggPrediction { diff --git a/image_analysis/bragg_prediction/BraggPredictionGPU.cu b/image_analysis/bragg_prediction/BraggPredictionGPU.cu index 8e172867..615eb794 100644 --- a/image_analysis/bragg_prediction/BraggPredictionGPU.cu +++ b/image_analysis/bragg_prediction/BraggPredictionGPU.cu @@ -127,13 +127,7 @@ namespace { out.d = 1.0f / sqrtf(recip_sq); out.dist_ewald = dist_ewald; out.rlp = 1.0f; - float partiality = 1.0f; - if (C.still_partiality && C.profile_radius_recipA > 0.0f) { - const float sig_bw = C.bandwidth_sigma * fabsf(recip_z); - const float sigma2 = C.profile_radius_recipA * C.profile_radius_recipA + sig_bw * sig_bw; - partiality = expf(-0.5f * dist_ewald * dist_ewald / sigma2); - } - out.partiality = partiality; + out.partiality = 1.0f; out.zeta = 1.0f; out.image_scale_corr = 1.0f; return true; @@ -164,9 +158,7 @@ namespace { float high_res_A, float ewald_dist_cutoff, char centering, - float bandwidth_sigma, - bool still_partiality, - float profile_radius_recipA) { + float bandwidth_sigma) { KernelConsts kc{}; auto geom = experiment.GetDiffractionGeometry(); kc.det_width_pxl = static_cast(experiment.GetXPixelsNum()); @@ -179,8 +171,6 @@ namespace { kc.one_over_wavelength = 1.0f / geom.GetWavelength_A(); kc.ewald_cutoff = ewald_dist_cutoff; kc.bandwidth_sigma = bandwidth_sigma; - kc.still_partiality = still_partiality; - kc.profile_radius_recipA = profile_radius_recipA; kc.Astar = lattice.Astar(); kc.Bstar = lattice.Bstar(); kc.Cstar = lattice.Cstar(); @@ -203,8 +193,7 @@ int BraggPredictionGPU::Calc(const DiffractionExperiment &experiment, const BraggPredictionSettings &settings) { // Build constants on host KernelConsts hK = BuildKernelConsts(experiment, lattice, settings.high_res_A, settings.ewald_dist_cutoff, - settings.centering, settings.bandwidth_sigma, - settings.still_partiality, settings.profile_radius_recipA); + settings.centering, settings.bandwidth_sigma); cudaMemcpyAsync(dK, &hK, sizeof(KernelConsts), cudaMemcpyHostToDevice, stream); cudaMemsetAsync(d_count, 0, sizeof(int), stream); diff --git a/image_analysis/bragg_prediction/BraggPredictionGPU.h b/image_analysis/bragg_prediction/BraggPredictionGPU.h index a7d2032d..ba0fc1f5 100644 --- a/image_analysis/bragg_prediction/BraggPredictionGPU.h +++ b/image_analysis/bragg_prediction/BraggPredictionGPU.h @@ -19,8 +19,6 @@ struct KernelConsts { float one_over_dmax_sq; float ewald_cutoff; float bandwidth_sigma; // relative Δλ/λ (sigma); 0 = monochromatic - bool still_partiality; // experimental stills excitation-error partiality (off => p = 1, fulls) - float profile_radius_recipA; // Gaussian sigma [1/A] for the stills partiality Coord Astar, Bstar, Cstar, S0; float rot[9]; char centering; diff --git a/image_analysis/indexing/AnalyzeIndexing.cpp b/image_analysis/indexing/AnalyzeIndexing.cpp index bc1661b8..af17244a 100644 --- a/image_analysis/indexing/AnalyzeIndexing.cpp +++ b/image_analysis/indexing/AnalyzeIndexing.cpp @@ -375,7 +375,10 @@ bool AnalyzeIndexing(DataMessage &message, int64_t indexing_lattice_count = 0; bool outcome = false; - const float min_frac = experiment.GetIndexingSettings().GetMinIndexedSpotFraction(); + // Minimum fraction of the in-resolution spots a candidate lattice must index to be accepted. + // Lowering it admits weaker/sparser crystals (more real ones on flooded XFEL frames, but also more + // spurious lattices that a downstream merge-consistency gate must remove). + constexpr float min_frac = 0.20f; if (nspots_indexed >= viable_cell_min_spots && nspots_indexed >= std::lround(min_frac * nspots_ref)) { auto uc = latt.GetUnitCell(); if (ok(uc.a) && ok(uc.b) && ok(uc.c) && ok(uc.alpha) && ok(uc.beta) && ok(uc.gamma)) { diff --git a/image_analysis/rotation_indexer/RotationIndexer.cpp b/image_analysis/rotation_indexer/RotationIndexer.cpp index 98bbecf8..20cf5603 100644 --- a/image_analysis/rotation_indexer/RotationIndexer.cpp +++ b/image_analysis/rotation_indexer/RotationIndexer.cpp @@ -253,9 +253,10 @@ void RotationIndexer::RunIndexing() { // Adopt the free triclinic cell only if it indexes CLEARLY more than the constrained cell - // a false promotion (a near-90 pseudo cell forced to ideal angles + a bogus centering) - // misplaces most reflections (CQ066 ratio ~0.1), whereas genuine higher symmetry (incl. - // R-centred) indexes comparably (ratio ~0.7). Preferring the constrained cell on a near-tie - // keeps the real symmetry/centering; the intensities settle the final space group. + // misplaces most reflections (measured indexed-fraction ratio ~0.1), whereas genuine higher + // symmetry (incl. R-centred) indexes comparably (ratio ~0.7). Preferring the constrained + // cell on a near-tie keeps the real symmetry/centering; the intensities settle the final + // space group. if (work[ci].has_tri) { Solved t = tri_f[ci].get(); if (t.ok && t.frac > 0.3f && frac < 0.5f * t.frac) { diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index fea7100b..52c31ed5 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -40,9 +40,7 @@ MergeOnTheFly::MergeOnTheFly(const DiffractionExperiment &x) high_resolution_limit(scaling_settings.GetHighResolutionLimit_A()), // A min-image-CC of 0 (the default) means "no limit": leave the optional // empty so the per-image CC cut is inactive. Otherwise a 0.0 threshold - // would silently drop every image with a non-positive per-image CC (which - // also wrongly zeroed N_obs in MergeStats, since it masks with cc_mask=true - // while the merge keeps all images). + // would silently drop every image with a non-positive per-image CC. image_cc_limit(scaling_settings.GetMinCCForImage() > 0.0 ? std::optional(scaling_settings.GetMinCCForImage()) : std::nullopt), @@ -64,10 +62,10 @@ bool MergeOnTheFly::IsMaskedRing(const Reflection &r) const { return ring >= 0 && ring < static_cast(masked_ice_rings.size()) && masked_ice_rings[ring]; } -void MergeOnTheFly::AddImage(const IntegrationOutcome &outcome, int64_t image_id, bool cc_mask) { +void MergeOnTheFly::AddImage(const IntegrationOutcome &outcome, int64_t image_id) { std::unique_lock ul(merged_mutex); - if (Mask(outcome, cc_mask)) + if (Mask(outcome)) return; const int half = HalfForImage(image_id); @@ -186,7 +184,7 @@ void MergeOnTheFly::RefineErrorModel(const std::vector &outc std::unordered_map> groups; for (const auto &outcome: outcomes) { - if (Mask(outcome, false)) + if (Mask(outcome)) continue; for (const auto &r: outcome.reflections) { if (generator.IsSystematicallyAbsent(r)) @@ -334,7 +332,7 @@ void MergeOnTheFly::RefineErrorModel(const std::vector &outc error_model_chi2 = chi2.empty() ? 0.0 : median(chi2) / CHI2_1_MEDIAN; } -bool MergeOnTheFly::Mask(const IntegrationOutcome &outcome, bool cc_mask) { +bool MergeOnTheFly::Mask(const IntegrationOutcome &outcome) { if (reference_cell) { auto cell = outcome.latt.GetUnitCell(); if (!cell.is_close(*reference_cell, @@ -343,7 +341,7 @@ bool MergeOnTheFly::Mask(const IntegrationOutcome &outcome, bool cc_mask) { return true; } - if (cc_mask && image_cc_limit) { + if (filter_by_image_cc && image_cc_limit) { if (!outcome.image_scale_cc || std::isnan(outcome.image_scale_cc.value()) || outcome.image_scale_cc.value() < image_cc_limit.value()) @@ -393,11 +391,10 @@ std::vector MergeOnTheFly::ExportReflections() { } std::vector MergeAll(const DiffractionExperiment &x, - const std::vector &integration_outcome, - bool mask) { + const std::vector &integration_outcome) { MergeOnTheFly merge(x); for (size_t i = 0; i < integration_outcome.size(); ++i) - merge.AddImage(integration_outcome[i], static_cast(i), mask); + merge.AddImage(integration_outcome[i], static_cast(i)); return merge.ExportReflections(); } @@ -413,6 +410,64 @@ struct ShellAccum { CorrelationCoefficient cc_ref; }; +std::pair ImageReferenceCC(const std::vector &reflections, + const std::map &reference, + const HKLKeyGenerator &generator, + std::optional d_min_limit, + double min_partiality) { + constexpr size_t MIN_REFLECTIONS = 20; + + double sum_x = 0.0; + double sum_y = 0.0; + double sum_x2 = 0.0; + double sum_y2 = 0.0; + double sum_xy = 0.0; + size_t n = 0; + + for (const auto &r: reflections) { + if (r.on_ice_ring) + continue; + if (!AcceptReflection(r, d_min_limit)) + continue; + if (r.partiality < min_partiality) + continue; + if (!std::isfinite(r.I) || !std::isfinite(r.image_scale_corr) || r.image_scale_corr <= 0.0f) + continue; + if (!std::isfinite(r.sigma) || r.sigma <= 0.0f) + continue; + + const auto it = reference.find(generator(r)); + if (it == reference.end()) + continue; + + const double image_i = static_cast(r.I) * static_cast(r.image_scale_corr); + const double ref_i = it->second; + + if (!std::isfinite(image_i) || !std::isfinite(ref_i)) + continue; + + sum_x += image_i; + sum_y += ref_i; + sum_x2 += image_i * image_i; + sum_y2 += ref_i * ref_i; + sum_xy += image_i * ref_i; + ++n; + } + + if (n < MIN_REFLECTIONS) + return {NAN, n}; + + const double nd = static_cast(n); + const double cov = sum_xy - sum_x * sum_y / nd; + const double var_x = sum_x2 - sum_x * sum_x / nd; + const double var_y = sum_y2 - sum_y * sum_y / nd; + + if (!(var_x > 0.0 && var_y > 0.0)) + return {NAN, n}; + + return {cov / std::sqrt(var_x * var_y), n}; +} + void CalcPossibleReflections(int space_group_number , const UnitCell &cell, double d_min, @@ -551,7 +606,7 @@ MergeStatistics MergeOnTheFly::MergeStats(const std::vector &m rmeas_obs.reserve(merged.size()); for (int i = 0; i < integration_outcome.size(); ++i) { - if (Mask(integration_outcome[i], true)) + if (Mask(integration_outcome[i])) continue; for (const auto &r: integration_outcome[i].reflections) { diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index cd19b7bd..a24509d4 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -100,6 +100,9 @@ class MergeOnTheFly { std::optional reference_cell; std::optional high_resolution_limit; std::optional image_cc_limit; + // Apply image_cc_limit in Mask(). One flag for the whole engine, not a per-call argument, so the + // merge, the error model and MergeStats can never disagree about which images are in. + bool filter_by_image_cc = false; double min_partiality = 0.02; // When set, ice-ring-flagged reflections are left out of this merge. Used for the P1 pass whose @@ -142,12 +145,13 @@ class MergeOnTheFly { std::unordered_map reject_median_I; size_t reject_count = 0; - bool Mask(const IntegrationOutcome &outcome, bool cc_mask); + bool Mask(const IntegrationOutcome &outcome); [[nodiscard]] bool IsMaskedRing(const Reflection &r) const; public: MergeOnTheFly(const DiffractionExperiment &x); MergeOnTheFly& ReferenceCell(const std::optional &cell); MergeOnTheFly& ExcludeIceRings(bool input) { exclude_ice_rings = input; return *this; } + MergeOnTheFly& FilterByImageCC(bool input) { filter_by_image_cc = input; return *this; } MergeOnTheFly& MaskIceRings(std::vector masked, float half_width_q) { masked_ice_rings = std::move(masked); mask_ice_half_width_q = half_width_q; return *this; } @@ -166,7 +170,7 @@ public: // image_id is the image's stable identity (its index in the outcomes vector). The CC1/2 half-set // is a deterministic hash of it, so the split is reproducible run-to-run and independent of the // order (or threading) of AddImage calls - not a draw from a shared RNG in call order. - void AddImage(const IntegrationOutcome& outcome, int64_t image_id, bool cc_mask = false); + void AddImage(const IntegrationOutcome& outcome, int64_t image_id); // d_min_override, when set, is the effective high-resolution limit for the shell table (used for // the automatic resolution cutoff computed by the caller); otherwise the manual @@ -180,5 +184,16 @@ public: }; std::vector MergeAll(const DiffractionExperiment &x, - const std::vector &reflections, - bool mask = false); + const std::vector &reflections); + +// Pearson CC between one image's corrected intensities (I * image_scale_corr) and a reference set of +// full intensities, over the reflections that would enter the merge (non-ice, within the resolution +// limit, partiality above the floor, finite). {NAN, n} when fewer than 20 reflections qualify. +// This is the per-image image_scale_cc: ScaleOnTheFly sets it, and StillsPartialityRefine recomputes it +// after refining the partiality model, so the reported CC always describes the corrections that will be +// merged - which matters because --min-image-cc drops images by it. +std::pair ImageReferenceCC(const std::vector &reflections, + const std::map &reference, + const HKLKeyGenerator &generator, + std::optional d_min_limit, + double min_partiality); diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 324a6c40..290bef99 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -169,7 +169,7 @@ RotationScaleMerge::RotationScaleMerge(const DiffractionExperiment &experiment, rfree_fraction = s.GetRfreeFraction(); scale_fulls = s.GetScaleFulls(); // Decay + absorption correction surfaces are one master toggle (on by default; both cross-validated, - // so a no-op when their systematic is absent). Decoupled from the stills-only -B / RefineB flag. + // so a no-op when their systematic is absent). refine_decay_b = s.GetCorrectionSurfaces(); absorption_iter = s.GetCorrectionSurfaces() ? s.GetAbsorptionIter() : 0; modulation_iter = s.GetCorrectionSurfaces() ? s.GetAbsorptionIter() : 0; @@ -1293,7 +1293,6 @@ void RotationScaleMerge::FinalizePerFrameScale(const std::vector &cc, co o.image_scale_cc_n.reset(); o.mosaicity_deg.reset(); } - o.image_scale_b_factor_Ang2.reset(); o.image_scale_wedge_deg.reset(); } } diff --git a/image_analysis/scale_merge/ScaleOnTheFly.cpp b/image_analysis/scale_merge/ScaleOnTheFly.cpp index 175b2aac..2c941c65 100644 --- a/image_analysis/scale_merge/ScaleOnTheFly.cpp +++ b/image_analysis/scale_merge/ScaleOnTheFly.cpp @@ -8,7 +8,6 @@ #include #include #include -#include namespace { // Robust loss scale (in sigma units) for the per-image scale fit: a few outlier reflections @@ -70,30 +69,6 @@ namespace { // The fixed-partiality residual for the Ceres path (used only when the B-factor is refined): the // stored partiality is a constant, so the model is G * partiality * exp(-B/(4 d^2)) * (1/rlp) * Itrue. - struct IntensityFixedResidual { - IntensityFixedResidual(const Reflection &r, double Itrue, double sigma) - : Iobs(r.I), - Itrue(Itrue), - weight(SafeInv(sigma, 1.0)), - lp(SafeInv(r.rlp, 1.0)), - b_resolution_coeff(-SafeInv(4.0 * r.d * r.d, 0.0)), - partiality(r.partiality) { - } - - template - bool operator()(const T *const G, const T *const B, T *residual) const { - const T B_term = ceres::exp(B[0] * T(b_resolution_coeff)); - residual[0] = (G[0] * T(partiality) * B_term * T(lp) * Itrue - T(Iobs)) * T(weight); - return true; - } - - const double Iobs; - const double Itrue; - const double weight; - const double lp; - const double b_resolution_coeff; - const double partiality; - }; } ScaleOnTheFly::ScaleOnTheFly(const DiffractionExperiment &x, const std::vector &ref) @@ -112,80 +87,23 @@ bool ScaleOnTheFly::Accept(const Reflection &r) const { return AcceptReflection(r, s.GetHighResolutionLimit_A()); } -std::pair ScaleOnTheFly::CalculateGlobalCC(const std::vector &reflections) const { - double sum_x = 0.0; - double sum_y = 0.0; - double sum_x2 = 0.0; - double sum_y2 = 0.0; - double sum_xy = 0.0; - size_t n = 0; - - for (const auto &r: reflections) { - if (r.on_ice_ring) - continue; - if (!AcceptReflection(r, s.GetHighResolutionLimit_A())) - continue; - if (r.partiality < s.GetMinPartiality()) - continue; - if (!std::isfinite(r.I) || !std::isfinite(r.image_scale_corr) || r.image_scale_corr <= 0.0f) - continue; - if (!std::isfinite(r.sigma) || r.sigma <= 0.0f) - continue; - - const HKLKey key = hkl_key_generator(r); - const auto it = reference_data.find(key); - if (it == reference_data.end()) - continue; - - const double image_i = static_cast(r.I) * static_cast(r.image_scale_corr); - const double ref_i = it->second; - - if (!std::isfinite(image_i) || !std::isfinite(ref_i)) - continue; - - sum_x += image_i; - sum_y += ref_i; - sum_x2 += image_i * image_i; - sum_y2 += ref_i * ref_i; - sum_xy += image_i * ref_i; - ++n; - } - - if (n < MIN_REFLECTIONS) - return {NAN, n}; - - const double nd = static_cast(n); - const double cov = sum_xy - sum_x * sum_y / nd; - const double var_x = sum_x2 - sum_x * sum_x / nd; - const double var_y = sum_y2 - sum_y * sum_y / nd; - - if (!(var_x > 0.0 && var_y > 0.0)) - return {NAN, n}; - - return {cov / std::sqrt(var_x * var_y), n}; -} - void ScaleOnTheFly::Scale(IntegrationOutcome &integration_outcome) const { if (integration_outcome.reflections.empty()) return; auto start = std::chrono::steady_clock::now(); - ScaleOnTheFlyResult result{ .B = 0.0, .G = 1.0 }; + ScaleOnTheFlyResult result{ .G = 1.0 }; auto clear_scale = [&]() { integration_outcome.image_scale_cc.reset(); integration_outcome.image_scale_cc_n.reset(); integration_outcome.image_scale_g.reset(); - integration_outcome.image_scale_b_factor_Ang2.reset(); }; - // With B fixed the fixed-partiality model G * coeff is linear in G, so the robust per-image scale is a - // 1-D M-estimate solved directly (IRLS) instead of a Ceres problem per image. Ceres is kept only when - // the B-factor (exp(-B/...)) is refined. - const bool linear_in_g = !s.GetRefineB(); - - if (linear_in_g) { + // The fixed-partiality model G * coeff is linear in G, so the robust per-image scale is a 1-D + // M-estimate solved directly (IRLS) rather than a Ceres problem per image. + { std::vector obs; obs.reserve(integration_outcome.reflections.size()); for (const auto &r: integration_outcome.reflections) { @@ -194,8 +112,7 @@ void ScaleOnTheFly::Scale(IntegrationOutcome &integration_outcome) const { const auto it = reference_data.find(hkl_key_generator(r)); if (it == reference_data.end()) continue; - const double B_term = std::exp(result.B * -SafeInv(4.0 * r.d * r.d, 0.0)); - const double coeff = r.partiality * B_term * SafeInv(r.rlp, 1.0) * it->second; + const double coeff = r.partiality * SafeInv(r.rlp, 1.0) * it->second; obs.push_back({coeff, static_cast(r.I), SafeInv(r.sigma, 1.0)}); } @@ -205,52 +122,18 @@ void ScaleOnTheFly::Scale(IntegrationOutcome &integration_outcome) const { } result.G = SolveScaleIRLS(obs, SCALE_ROBUST_K); - } else { - ceres::Problem problem; - - size_t n_reflections = 0; - for (const auto &r: integration_outcome.reflections) { - if (!Accept(r)) - continue; - - const HKLKey key = hkl_key_generator(r); - if (!reference_data.contains(key)) - continue; - - ++n_reflections; - - auto *cost = new ceres::AutoDiffCostFunction( - new IntensityFixedResidual(r, reference_data.at(key), r.sigma)); - problem.AddResidualBlock(cost, new ceres::CauchyLoss(SCALE_ROBUST_K), &result.G, &result.B); - } - - if (n_reflections < MIN_REFLECTIONS) { - clear_scale(); - return; - } - - problem.SetParameterLowerBound(&result.G, 0, 0.0); - problem.SetParameterLowerBound(&result.B, 0, s.GetMinB()); - problem.SetParameterUpperBound(&result.B, 0, s.GetMaxB()); - - ceres::Solver::Options options; - options.linear_solver_type = ceres::DENSE_QR; - options.minimizer_progress_to_stdout = false; - options.num_threads = 1; - - ceres::Solver::Summary summary; - ceres::Solve(options, &problem, &summary); } for (auto &r: integration_outcome.reflections) { - const double B_term = exp(result.B * -SafeInv(4.0 * r.d * r.d, 0.0)); - const double denom = B_term * r.partiality * result.G; + const double denom = r.partiality * result.G; r.image_scale_corr = (std::isfinite(r.rlp) && std::isfinite(denom) && denom > 0.0) ? static_cast(r.rlp / denom) : NAN; } - const auto [cc, cc_n] = CalculateGlobalCC(integration_outcome.reflections); + const auto [cc, cc_n] = ImageReferenceCC(integration_outcome.reflections, reference_data, + hkl_key_generator, s.GetHighResolutionLimit_A(), + s.GetMinPartiality()); result.cc = cc; result.cc_n = cc_n; @@ -261,11 +144,6 @@ void ScaleOnTheFly::Scale(IntegrationOutcome &integration_outcome) const { integration_outcome.image_scale_cc_n = cc_n; integration_outcome.image_scale_g = result.G; integration_outcome.image_scale_wedge_deg.reset(); - - if (s.GetRefineB()) - integration_outcome.image_scale_b_factor_Ang2 = result.B; - else - integration_outcome.image_scale_b_factor_Ang2.reset(); } void ScaleOnTheFly::Scale(std::vector &integration, size_t nthreads) const { diff --git a/image_analysis/scale_merge/ScaleOnTheFly.h b/image_analysis/scale_merge/ScaleOnTheFly.h index 67bd5623..bc8ae639 100644 --- a/image_analysis/scale_merge/ScaleOnTheFly.h +++ b/image_analysis/scale_merge/ScaleOnTheFly.h @@ -13,7 +13,6 @@ struct ScaleOnTheFlyResult { - double B = 0; double G = 1.0; double cc = NAN; size_t cc_n = 0; @@ -23,9 +22,8 @@ struct ScaleOnTheFlyResult { // Per-image reference scaling with the FIXED partiality model: each reflection's stored partiality is // used as-is (it is 1 for stills and the zeta/erf rocking-curve value already set at prediction for -// rotation data). No partiality recompute, no mosaicity or wedge refinement. The B-factor may be refined -// (Ceres); otherwise the fit is linear in G (a robust 1-D IRLS). Rotation self-scaling/merging offline -// uses the dedicated RotationScaleMerge path instead. +// rotation data). No partiality recompute, no mosaicity or wedge refinement: the fit is linear in G, a +// robust 1-D IRLS. Rotation self-scaling/merging offline uses the dedicated RotationScaleMerge path. class ScaleOnTheFly { constexpr static size_t MIN_REFLECTIONS = 20; @@ -34,7 +32,6 @@ class ScaleOnTheFly { std::map reference_data; bool Accept(const Reflection &r) const; - [[nodiscard]] std::pair CalculateGlobalCC(const std::vector &reflections) const; public: ScaleOnTheFly(const DiffractionExperiment &x, const std::vector &ref); diff --git a/image_analysis/scale_merge/ScalingResult.cpp b/image_analysis/scale_merge/ScalingResult.cpp index 83874371..1a02d09e 100644 --- a/image_analysis/scale_merge/ScalingResult.cpp +++ b/image_analysis/scale_merge/ScalingResult.cpp @@ -10,7 +10,6 @@ ScalingResult::ScalingResult(size_t n) : image_scale_g(n, NAN), mosaicity_deg(n, NAN), - image_bfactor_Ang2(n, NAN), rotation_wedge_deg(n, NAN), image_cc(n, NAN), image_cc_n(n, 0) {} @@ -18,14 +17,12 @@ ScalingResult::ScalingResult(size_t n) ScalingResult::ScalingResult(const std::vector &v) : image_scale_g(v.size(), NAN), mosaicity_deg(v.size(), NAN), - image_bfactor_Ang2(v.size(), NAN), rotation_wedge_deg(v.size(), NAN), image_cc(v.size(), NAN), image_cc_n(v.size(), 0) { for (int i = 0; i < v.size(); i++) { image_scale_g[i] = v[i].image_scale_g.value_or(NAN); mosaicity_deg[i] = v[i].mosaicity_deg.value_or(NAN); - image_bfactor_Ang2[i] = v[i].image_scale_b_factor_Ang2.value_or(NAN); rotation_wedge_deg[i] = v[i].image_scale_wedge_deg.value_or(NAN); image_cc[i] = v[i].image_scale_cc.value_or(NAN); image_cc_n[i] = v[i].image_scale_cc_n.value_or(0); @@ -41,12 +38,11 @@ void ScalingResult::SaveToFile(const std::string &filename) { } // Header so the columns are self-describing (lines starting with '#' are comments for gnuplot/numpy). - img_file << "# image_number scale_G bfactor_Ang2 mosaicity_deg wedge_deg cc_to_merge cc_n\n"; + img_file << "# image_number scale_G mosaicity_deg wedge_deg cc_to_merge cc_n\n"; for (size_t i = 0; i < image_scale_g.size(); ++i) { img_file << i << " " << image_scale_g[i] - << " " << image_bfactor_Ang2[i] << " " << mosaicity_deg[i] << " " << rotation_wedge_deg[i] << " " << image_cc[i] diff --git a/image_analysis/scale_merge/ScalingResult.h b/image_analysis/scale_merge/ScalingResult.h index 7003b5a8..717d1915 100644 --- a/image_analysis/scale_merge/ScalingResult.h +++ b/image_analysis/scale_merge/ScalingResult.h @@ -10,7 +10,6 @@ struct ScalingResult { std::vector image_scale_g; std::vector mosaicity_deg; - std::vector image_bfactor_Ang2; std::vector rotation_wedge_deg; std::vector image_cc; std::vector image_cc_n; diff --git a/image_analysis/scale_merge/StillsPartialityRefine.cpp b/image_analysis/scale_merge/StillsPartialityRefine.cpp index 03f60189..9387fdff 100644 --- a/image_analysis/scale_merge/StillsPartialityRefine.cpp +++ b/image_analysis/scale_merge/StillsPartialityRefine.cpp @@ -39,7 +39,11 @@ namespace { // Analytic partiality for a reflection whose base reciprocal vector is q, tilted by (psi_x, psi_y). // Mirrors BraggPrediction: dist_ewald = |S| - 1/lambda with S = q_rot + S0, and - // p = exp(-dist_ewald^2 / 2 sigma^2), sigma^2 = (gamma0 + gamma_e*d*)^2 + (bw*|q_z|)^2. + // p = exp(-dist_ewald^2 / 2 sigma^2), sigma^2 = gamma0^2 + (gamma_e*d*)^2 + (bw*|q_z|)^2 - three + // independent broadenings added in quadrature: the reciprocal-lattice point's own radius (gamma0, + // ~1/domain size, resolution-INdependent), the mosaic/divergence spread (gamma_e*d*, proportional to + // d*) and the bandwidth smear along the beam. In practice the fit below returns gamma_e ~ 0 and the + // width is essentially gamma0 - see there. double ComputeP(double qx, double qy, double qz, double psi_x, double psi_y, double s0x, double s0y, double s0z, double inv_lambda, @@ -51,9 +55,9 @@ namespace { const double Sx = qr[0] + s0x, Sy = qr[1] + s0y, Sz = qr[2] + s0z; const double de = std::sqrt(Sx * Sx + Sy * Sy + Sz * Sz) - inv_lambda; const double dstar = std::sqrt(qr[0] * qr[0] + qr[1] * qr[1] + qr[2] * qr[2]); - const double sig = gamma0 + gamma_e * dstar; + const double sig_ang = gamma_e * dstar; const double sbw = bw * std::fabs(qr[2]); - const double sig2 = sig * sig + sbw * sbw; + const double sig2 = gamma0 * gamma0 + sig_ang * sig_ang + sbw * sbw; if (!(sig2 > 0.0)) return 1.0; return std::exp(-0.5 * de * de / sig2); @@ -113,9 +117,9 @@ namespace { const T Sx = qr[0] + T(s0x), Sy = qr[1] + T(s0y), Sz = qr[2] + T(s0z); const T de = ceres::sqrt(Sx * Sx + Sy * Sy + Sz * Sz) - T(inv_lambda); const T dstar = ceres::sqrt(qr[0] * qr[0] + qr[1] * qr[1] + qr[2] * qr[2]); - const T sig = T(gamma0) + T(gamma_e) * dstar; + const T sig_ang = T(gamma_e) * dstar; const T sbw = T(bw) * ceres::abs(qr[2]); - const T sig2 = sig * sig + sbw * sbw; + const T sig2 = T(gamma0) * T(gamma0) + sig_ang * sig_ang + sbw * sbw; const T p = ceres::exp(T(-0.5) * de * de / sig2); residual[0] = T(weight) * (T(G) * p * T(lp) * T(Iref) - T(Iobs)); return true; @@ -136,13 +140,10 @@ namespace { } StillsPartialityRefine::StillsPartialityRefine(const DiffractionExperiment &x) - : StillsPartialityRefine(x, Settings{}) {} - -StillsPartialityRefine::StillsPartialityRefine(const DiffractionExperiment &x, Settings settings) : experiment_(x), - settings_(settings), hkl_key_generator_(x.GetScalingSettings().GetMergeFriedel(), x.GetSpaceGroupNumber().value_or(1)), d_min_limit_(x.GetScalingSettings().GetHighResolutionLimit_A()), + min_partiality_(x.GetScalingSettings().GetMinPartiality()), bandwidth_sigma_(x.GetBandwidthFWHM().value_or(0.0f) / 2.3548f) {} double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome, @@ -156,7 +157,6 @@ double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome, const Coord S0 = outcome.geom.GetScatteringVector(); const double inv_lambda = 1.0 / outcome.geom.GetWavelength_A(); const double bw = bandwidth_sigma_; - const double gamma0 = 0.0; // width is purely angular: sigma(d*) = gamma_e * d* (set per crystal below) auto base_q = [&](const Reflection &r) { return Astar * static_cast(r.h) + Bstar * static_cast(r.k) @@ -166,7 +166,9 @@ double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome, // Collect the reflections that constrain the fit (accepted, non-ice, finite, present in the reference). std::vector obs; obs.reserve(outcome.reflections.size()); - double sum_ang2 = 0.0; // RMS angular excitation error (dist_ewald / d*) -> per-crystal mosaic width + // Moments of the excitation error against resolution: de^2 ~ gamma0^2 + gamma_e^2 * d*^2, fitted per + // crystal by ordinary least squares on (d*^2, de^2). Both components come out of the data. + double m_n = 0.0, m_x = 0.0, m_xx = 0.0, m_y = 0.0, m_xy = 0.0; size_t n_de = 0; for (const Reflection &r: outcome.reflections) { if (r.on_ice_ring || !AcceptReflection(r, d_min_limit_)) @@ -186,15 +188,20 @@ double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome, .weight = SafeInv(r.sigma, 1.0), }); - // Angular excitation error delta_psi = dist_ewald / d* at the stored orientation (psi = 0). Using - // the ANGULAR distance (not the linear reciprocal-space distance) makes the partiality width - // resolution-clean: a fixed mosaic angle smears high-resolution rlps more in reciprocal space, so a - // constant linear width computes p too small at high resolution and over-divides those shells. + // Excitation error at the stored orientation (psi = 0), collected as the moments of de^2 against + // d*^2, so BOTH width components are fitted rather than one being forced to zero. Forcing the + // width to be purely angular (gamma0 = 0) pins it to the high-resolution edge - it is fitted over + // a d*^2-dense population - and it then collapses at low d*, giving p ~ 0 for reflections that + // were plainly recorded, which inflated the merged low-resolution intensity scale ~3.6x. const double dstar = VecLen(q.x, q.y, q.z); const double de0 = VecLen(q.x + S0.x, q.y + S0.y, q.z + S0.z) - inv_lambda; if (dstar > 1e-9) { - const double dpsi = de0 / dstar; - sum_ang2 += dpsi * dpsi; + const double x = dstar * dstar, y = de0 * de0; + m_n += 1.0; + m_x += x; + m_xx += x * x; + m_y += y; + m_xy += x * y; ++n_de; } } @@ -202,12 +209,28 @@ double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome, if (obs.size() < MIN_FIT_REFLECTIONS || n_de == 0) return 0.0; - // Per-crystal angular mosaic width from the RMS angular excitation error. sigma(d*) = gamma_e * d* - // (gamma0 = 0), i.e. p = exp(-0.5 (delta_psi / gamma_e)^2) is a Gaussian in the angular distance from - // the Ewald sphere - the physical mosaic/divergence model, independent of resolution. A positive - // settings_.gamma_e overrides the per-crystal estimate with a shared (pooled) width. - const double gamma_e_ang = std::max(std::sqrt(sum_ang2 / static_cast(n_de)), 1e-9); - const double gamma_e = settings_.gamma_e > 0.0 ? settings_.gamma_e : gamma_e_ang; + // Solve the 2x2 normal equations for de^2 = A + B d*^2. A degenerate spread in d* (all reflections in + // one shell) leaves B undetermined, so fall back to the pure angular width there; a negative fitted + // component is unphysical and is clamped to zero, which reduces to the previous model. + // + // Measured outcome, worth knowing before touching this: the fit does NOT split the width between the + // two terms - it returns gamma0 ~ 4e-4 1/A and gamma_e ~ 0 (their cross-over sits at d = 0.3 A, far + // outside any measured range), i.e. a width constant in the LINEAR Ewald distance. That is + // structural, not a fluke: prediction accepts reflections on a fixed linear |dist_ewald| cutoff, so + // the accepted population's de^2 is flat in d*^2 by construction and the slope is genuinely ~0. The + // truncated population cannot constrain an angular term; the resolution-independent one is what the + // data actually support. + const double det = m_n * m_xx - m_x * m_x; + double A = 0.0, B = 0.0; + if (std::fabs(det) > 1e-30) { + A = (m_xx * m_y - m_x * m_xy) / det; + B = (m_n * m_xy - m_x * m_y) / det; + } else { + B = m_x > 0.0 ? m_y / m_x : 0.0; + } + const double gamma0 = std::sqrt(std::max(0.0, A)); + const double gamma_e_fit = std::max(std::sqrt(std::max(0.0, B)), 1e-9); + const double gamma_e = settings_.gamma_e > 0.0 ? settings_.gamma_e : gamma_e_fit; double psi[2] = {0.0, 0.0}; double G = 1.0; @@ -262,19 +285,28 @@ double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome, } // Write the refined partiality + scale correction onto every reflection of the crystal (not only the - // fit subset), so the merge sees a consistent model. image_scale_corr = rlp / (partiality * G). + // fit subset), so the merge sees a consistent model. image_scale_corr = rlp / (partiality * G), the + // same composition ScaleOnTheFly writes. for (auto &r: outcome.reflections) { const Coord q = base_q(r); const double p = ComputeP(q.x, q.y, q.z, psi[0], psi[1], S0.x, S0.y, S0.z, inv_lambda, gamma0, gamma_e, bw); r.partiality = static_cast(p); const double denom = p * G; - r.image_scale_corr = (std::isfinite(r.rlp) && denom > 0.0) + r.image_scale_corr = (std::isfinite(r.rlp) && std::isfinite(denom) && denom > 0.0) ? static_cast(r.rlp / denom) : NAN; } outcome.image_scale_g = static_cast(G); + // The corrections just changed, so the CC that ScaleOnTheFly measured no longer describes them. + // Refresh it here: it is reported per image and --min-image-cc drops images by it, so it has to be + // the CC of the data that is actually merged. + const auto [cc, cc_n] = ImageReferenceCC(outcome.reflections, reference, hkl_key_generator_, + d_min_limit_, min_partiality_); + outcome.image_scale_cc = cc; + outcome.image_scale_cc_n = cc_n; + const double tilt_deg = std::sqrt(psi[0] * psi[0] + psi[1] * psi[1]) * kRadToDeg; return tilt_deg; } @@ -288,7 +320,7 @@ double StillsPartialityRefine::Run(std::vector &outcomes, si for (int outer = 0; outer < settings_.outer_iterations; ++outer) { // Reference full intensities from the current corrections. - const std::vector merged = MergeAll(experiment_, outcomes, false); + const std::vector merged = MergeAll(experiment_, outcomes); std::map reference; for (const auto &m: merged) reference[hkl_key_generator_(m)] = m.I; diff --git a/image_analysis/scale_merge/StillsPartialityRefine.h b/image_analysis/scale_merge/StillsPartialityRefine.h index 108f7cac..4db99990 100644 --- a/image_analysis/scale_merge/StillsPartialityRefine.h +++ b/image_analysis/scale_merge/StillsPartialityRefine.h @@ -9,7 +9,7 @@ #include "../IntegrationOutcome.h" #include "HKLKey.h" -// Experimental physical partiality post-refinement for STILLS (env JFJOCH_STILL_POSTREFINE). +// Physical partiality post-refinement for STILLS (on by default; rugnux --simple-stills opts out). // // The default stills partiality is a frozen scalar-sigma Gaussian (or p == 1): p is set once at // prediction and never optimised, and any attempt to free a per-image sigma jointly with the per-image @@ -41,16 +41,16 @@ public: }; explicit StillsPartialityRefine(const DiffractionExperiment &x); - StillsPartialityRefine(const DiffractionExperiment &x, Settings settings); // Refine all crystals in place. Returns the mean |dpsi| applied (degrees), for diagnostics. double Run(std::vector &outcomes, size_t nthreads = 0) const; private: const DiffractionExperiment experiment_; - const Settings settings_; + const Settings settings_{}; const HKLKeyGenerator hkl_key_generator_; const std::optional d_min_limit_; + const double min_partiality_; const float bandwidth_sigma_; // Refine one crystal against the reference map; returns |dpsi| in degrees (0 if skipped). diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp index a722a67c..8ccdcf06 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp +++ b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp @@ -4,7 +4,6 @@ #include #include #include -#include #include "AdaptiveSpotFinderCPU.h" #include "AdaptiveThreshold.h" @@ -59,9 +58,8 @@ void AdaptiveSpotFinderCPU::AccumulateRings(const ImagePreprocessorBuffer &image } } -std::vector AdaptiveSpotFinderCPU::Run(const ImagePreprocessorBuffer &image, - const SpotFindingSettings &settings, - const std::vector &res_mask) { +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; @@ -79,8 +77,11 @@ std::vector AdaptiveSpotFinderCPU::Run(const ImagePreprocessorB g_sum += ring_sum[b]; g_sum2 += ring_sum2[b]; } - if (n_total == 0) - return {}; + if (n_total == 0) { + // Nothing valid to threshold against: leave no strong pixels for ExtractSpots to build on. + std::fill(output_buffer.begin(), output_buffer.end(), 0); + return; + } const double E = std::max(1.0f, settings.false_pixels_per_frame); double p = E / static_cast(n_total); @@ -122,7 +123,4 @@ std::vector AdaptiveSpotFinderCPU::Run(const ImagePreprocessorB } if (npix % 32 != 0) output_buffer[OutputSize() - 1] = out.to_ulong(); - - // --- Stage D: connected components + resolution mask + min/max-pix (shared with classic path) --- - return ExtractSpots(image, settings, res_mask); } diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h index d75e937b..bb4aa220 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h +++ b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h @@ -39,7 +39,5 @@ class AdaptiveSpotFinderCPU : public ImageSpotFinder { public: explicit AdaptiveSpotFinderCPU(const AzimuthalIntegrationMapping &mapping); - std::vector Run(const ImagePreprocessorBuffer &image, - const SpotFindingSettings &settings, - const std::vector &res_mask) override; + void Detect(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) override; }; diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu index 4403bde3..3a0b4d5d 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu +++ b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu @@ -257,12 +257,11 @@ void AdaptiveSpotFinderGPU::ComputeThresholds(const SpotFindingSettings &setting } } -std::vector AdaptiveSpotFinderGPU::Run(const ImagePreprocessorBuffer &image, - const SpotFindingSettings &settings, - const std::vector &res_mask) { +void AdaptiveSpotFinderGPU::Detect(const ImagePreprocessorBuffer &image, + const SpotFindingSettings &settings) { if (image.size() != npix) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - "AdaptiveSpotFinderGPU::Run: mismatch in pixel size"); + "AdaptiveSpotFinderGPU::Detect: mismatch in pixel size"); // --- Stage A: robust per-ring background (one plain pass + two sigma-clip passes) --- cuda_err(cudaMemsetAsync(gpu_sum, 0, sizeof(float) * nbins, *stream)); @@ -303,8 +302,11 @@ std::vector AdaptiveSpotFinderGPU::Run(const ImagePreprocessorB last_profile.Clear(mapping); last_profile.Add(prof_sum, prof_sum2, prof_count); - if (host_thr.empty()) - return {}; + if (host_thr.empty()) { + // Nothing valid to threshold against: leave no strong pixels for ExtractSpots to build on. + std::fill(output_buffer.begin(), output_buffer.end(), 0); + return; + } // --- Stage C: flag strong pixels into the bit buffer (value >= ring threshold) --- cuda_err(cudaMemcpyAsync(gpu_thr, host_thr.data(), sizeof(float) * nbins, cudaMemcpyHostToDevice, *stream)); @@ -313,7 +315,4 @@ std::vector AdaptiveSpotFinderGPU::Run(const ImagePreprocessorB image.getGPUBuffer(), gpu_pixel_to_bin, gpu_thr, gpu_strong, npix, nbins); cuda_err(cudaMemcpyAsync(output_buffer.data(), gpu_strong, OutputByteSize(), cudaMemcpyDeviceToHost, *stream)); cuda_err(cudaStreamSynchronize(*stream)); - - // --- Stage D: connected components + resolution mask + min/max-pix (shared host path) --- - return ExtractSpots(image, settings, res_mask); } diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h index e115ca8f..78740f2a 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h +++ b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h @@ -93,11 +93,9 @@ public: AdaptiveSpotFinderGPU(const AdaptiveSpotFinderGPU &) = delete; AdaptiveSpotFinderGPU &operator=(const AdaptiveSpotFinderGPU &) = delete; - std::vector Run(const ImagePreprocessorBuffer &image, - const SpotFindingSettings &settings, - const std::vector &res_mask) override; + void Detect(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) override; - // The azimuthal profile computed as a byproduct of the last Run() - lets this engine replace the + // The azimuthal profile computed as a byproduct of the last Detect() - lets this engine replace the // separate azint pass in the analysis pipeline. [[nodiscard]] const AzimuthalIntegrationProfile &GetProfile() const { return last_profile; } }; diff --git a/image_analysis/spot_finding/DetModuleSpotFinder_cpu.h b/image_analysis/spot_finding/DetModuleSpotFinder_cpu.h index 5506387c..e0edbf2a 100644 --- a/image_analysis/spot_finding/DetModuleSpotFinder_cpu.h +++ b/image_analysis/spot_finding/DetModuleSpotFinder_cpu.h @@ -39,7 +39,7 @@ void FindSpots(DeviceOutput &output, || (col == 767) || (col == 768)) bad_pixel = 1; - if ((d_array[coord] < settings.high_resolution_limit) + if ((d_array[coord] < settings.high_resolution_limit.value_or(0.0f)) || (d_array[coord] > settings.low_resolution_limit)) bad_pixel = 1; diff --git a/image_analysis/spot_finding/ImageSpotFinder.cpp b/image_analysis/spot_finding/ImageSpotFinder.cpp index 895b8dba..9d0e2939 100644 --- a/image_analysis/spot_finding/ImageSpotFinder.cpp +++ b/image_analysis/spot_finding/ImageSpotFinder.cpp @@ -20,6 +20,13 @@ size_t ImageSpotFinder::OutputByteSize() const { return OutputSize() * sizeof(uint32_t); } +std::vector ImageSpotFinder::Run(const ImagePreprocessorBuffer &image, + const SpotFindingSettings &settings, + const std::vector &res_mask) { + Detect(image, settings); + return ExtractSpots(image, settings, res_mask); +} + std::vector ImageSpotFinder::ExtractSpots(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings, const std::vector &res_mask) { diff --git a/image_analysis/spot_finding/ImageSpotFinder.h b/image_analysis/spot_finding/ImageSpotFinder.h index 11691aee..b4b938a3 100644 --- a/image_analysis/spot_finding/ImageSpotFinder.h +++ b/image_analysis/spot_finding/ImageSpotFinder.h @@ -18,13 +18,19 @@ protected: ImageSpotFinder(int32_t width, int32_t height); size_t OutputSize() const; size_t OutputByteSize() const; - std::vector ExtractSpots(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings, const std::vector &res_mask); public: constexpr static int32_t MIN_VALID_PIXELS = 100; constexpr static int NBX = 15; virtual ~ImageSpotFinder() = default; - virtual std::vector Run(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings, const std::vector &res_mask) = 0; + + // Detect flags the image's strong pixels into the internal bit buffer - the expensive step (local + // box or per-ring background over every pixel). ExtractSpots then builds the spots from those + // pixels; min/max-pix and the resolution mask enter only there, so several min-pix values can be + // tried on ONE detection pass (MXAnalysisWithoutFPGA does that when min-pix is chosen per image). + virtual void Detect(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) = 0; + std::vector ExtractSpots(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings, const std::vector &res_mask); + std::vector Run(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings, const std::vector &res_mask); }; diff --git a/image_analysis/spot_finding/ImageSpotFinderCPU.cpp b/image_analysis/spot_finding/ImageSpotFinderCPU.cpp index 97188897..9b203baf 100644 --- a/image_analysis/spot_finding/ImageSpotFinderCPU.cpp +++ b/image_analysis/spot_finding/ImageSpotFinderCPU.cpp @@ -9,9 +9,8 @@ ImageSpotFinderCPU::ImageSpotFinderCPU(int32_t in_width, int32_t in_height) : ImageSpotFinder(in_width, in_height) {} -std::vector ImageSpotFinderCPU::Run(const ImagePreprocessorBuffer &image, - const SpotFindingSettings &settings, - const std::vector &res_mask) { +void ImageSpotFinderCPU::Detect(const ImagePreprocessorBuffer &image, + const SpotFindingSettings &settings) { for (int i = 0; i < OutputSize(); i++) output_buffer[i] = 0; @@ -131,6 +130,4 @@ std::vector ImageSpotFinderCPU::Run(const ImagePreprocessorBuff if (height * width % 32 != 0) output_buffer[OutputSize() - 1] = out.to_ulong(); - - return ExtractSpots(image, settings, res_mask); } diff --git a/image_analysis/spot_finding/ImageSpotFinderCPU.h b/image_analysis/spot_finding/ImageSpotFinderCPU.h index 29ae00f4..2f8c28c3 100644 --- a/image_analysis/spot_finding/ImageSpotFinderCPU.h +++ b/image_analysis/spot_finding/ImageSpotFinderCPU.h @@ -20,7 +20,7 @@ class ImageSpotFinderCPU : public ImageSpotFinder { public: ImageSpotFinderCPU(int32_t width, int32_t height); - std::vector Run(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings, const std::vector &res_mask); + void Detect(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) override; }; diff --git a/image_analysis/spot_finding/ImageSpotFinderGPU.cu b/image_analysis/spot_finding/ImageSpotFinderGPU.cu index f1afc258..2a1afc76 100644 --- a/image_analysis/spot_finding/ImageSpotFinderGPU.cu +++ b/image_analysis/spot_finding/ImageSpotFinderGPU.cu @@ -237,7 +237,7 @@ ImageSpotFinderGPU::ImageSpotFinderGPU(int32_t in_width, int32_t in_height, gpu_out_1 = CudaDevicePtr(OutputSize()); } -std::vector ImageSpotFinderGPU::Run(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings, const std::vector &res_mask) { +void ImageSpotFinderGPU::Detect(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) { spot_parameters spot_params{}; spot_params.height = height; spot_params.width = width; @@ -273,6 +273,4 @@ std::vector ImageSpotFinderGPU::Run(const ImagePreprocessorBuff cuda_err(cudaMemcpyAsync(output_buffer.data(), gpu_out_1, OutputSize() * sizeof(uint32_t), cudaMemcpyDeviceToHost, *stream)); cuda_err(cudaStreamSynchronize(*stream)); - - return ExtractSpots(image, settings, res_mask); } diff --git a/image_analysis/spot_finding/ImageSpotFinderGPU.h b/image_analysis/spot_finding/ImageSpotFinderGPU.h index c18afbe2..0dab3a74 100644 --- a/image_analysis/spot_finding/ImageSpotFinderGPU.h +++ b/image_analysis/spot_finding/ImageSpotFinderGPU.h @@ -23,7 +23,7 @@ public: ImageSpotFinderGPU(int32_t width, int32_t height, std::shared_ptr stream); ~ImageSpotFinderGPU() override = default; - std::vector Run(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings, const std::vector &res_mask) override; + void Detect(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) override; }; diff --git a/image_analysis/spot_finding/SpotFindingSettings.h b/image_analysis/spot_finding/SpotFindingSettings.h index 96a4c9af..76cff2b0 100644 --- a/image_analysis/spot_finding/SpotFindingSettings.h +++ b/image_analysis/spot_finding/SpotFindingSettings.h @@ -16,7 +16,10 @@ struct SpotFindingSettings { // 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; + // 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). + std::optional high_resolution_limit; float low_resolution_limit = 50.0; float cutoff_spot_count_low_res = 5.0; std::optional high_res_gap_Q_recipA = 1.5; // 0.25 * 2 * pi diff --git a/image_analysis/spot_finding/SpotUtils.cpp b/image_analysis/spot_finding/SpotUtils.cpp index b682649b..638e0c7c 100644 --- a/image_analysis/spot_finding/SpotUtils.cpp +++ b/image_analysis/spot_finding/SpotUtils.cpp @@ -147,7 +147,8 @@ void SpotAnalyze(const DiffractionExperiment &experiment, CountSpots(output, spots_out, spot_finding_settings.cutoff_spot_count_low_res); - GenerateSpotPlot(output, spots_out, spot_finding_settings.high_resolution_limit); + GenerateSpotPlot(output, spots_out, + spot_finding_settings.high_resolution_limit.value_or(experiment.GetDetectorMaxResolution_A())); output.resolution_estimate = GetResolution(spots_out); diff --git a/reader/HDF5MetadataSource.cpp b/reader/HDF5MetadataSource.cpp index 4a0cb2f3..75eac328 100644 --- a/reader/HDF5MetadataSource.cpp +++ b/reader/HDF5MetadataSource.cpp @@ -372,7 +372,6 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen dataset->b_factor = master_file->ReadOptVector("/entry/MX/bFactor"); dataset->image_scale_factor = master_file->ReadOptVector("/entry/MX/imageScaleFactor"); dataset->image_scale_cc = master_file->ReadOptVector("/entry/MX/imageScaleCC"); - dataset->image_scale_b = master_file->ReadOptVector("/entry/MX/imageScaleBFactor"); dataset->integrated_reflections = master_file->ReadOptVector("/entry/MX/integratedReflections"); } if (master_file->Exists("/entry/image")) @@ -951,8 +950,6 @@ void HDF5MetadataSource::FillPerImage(DataMessage &message, int64_t requested_im message.mosaicity_deg = dataset->mosaicity_deg[image_number]; if (dataset->b_factor.size() > image_number) message.b_factor = dataset->b_factor[image_number]; - if (dataset->image_scale_b.size() > image_number) - message.image_scale_b_factor = dataset->image_scale_b[image_number]; if (dataset->image_scale_factor.size() > image_number) message.image_scale_factor = dataset->image_scale_factor[image_number]; if (dataset->image_scale_cc.size() > image_number) diff --git a/reader/JFJochReaderDataset.h b/reader/JFJochReaderDataset.h index 5dcae787..c38ea2c1 100644 --- a/reader/JFJochReaderDataset.h +++ b/reader/JFJochReaderDataset.h @@ -49,7 +49,6 @@ struct JFJochReaderDataset { std::vector integrated_reflections; std::vector image_scale_factor; std::vector image_scale_cc; - std::vector image_scale_b; std::vector max_value; // Maps this dataset's image index -> the original/collected image number it came from. diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index b075daaf..f51c5e14 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -790,6 +790,13 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b result.cancelled = cancelled_; result.images_processed = finished_count.load(); + + // Every image failing is a total failure, not a run that produced nothing: it used to be reported + // only as per-image log lines while the process still exited 0 with no output file. + if (!cancelled_ && images_to_process > 0 && result.images_processed == 0) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "No image could be analyzed - see the per-image errors above"); + result.mean_processing_time = plots.GetMeanProcessingTime(); result.indexing_rate = plots.GetIndexingRate(); @@ -938,12 +945,11 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b const bool is_rotation = experiment_.IsRotationIndexing(); // rotation indexing -> rotation scaling/merge std::optional rsm; if (is_rotation) { - if (rot_ss.GetRefineB() - || experiment_.GetRefineRotationWedgeInScaling() + if (experiment_.GetRefineRotationWedgeInScaling() || rot_ss.GetRotationWedgeForScaling().has_value()) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Rotation scaling/merging (RotationScaleMerge) does not support " - "B-factor refinement or wedge refinement"); + "wedge refinement"); // A reference MTZ is allowed for rotation: it fixes the space group / cell (on the CLI) and // resolves the indexing ambiguity (below), but is NOT used to scale - the rotation merge stays // self-consistent. @@ -971,7 +977,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // degrade weak stills (measured CC1/2 collapse at the default 3 iters). { phase("Scaling images (" + label + ")"); - auto merge_result = MergeAll(experiment_, indexer->GetIntegrationOutcome(), false); + auto merge_result = MergeAll(experiment_, indexer->GetIntegrationOutcome()); indexer->ScaleAllImages(merge_result); } // Physical partiality post-refinement (default on; --simple-stills disables): refine a per-crystal @@ -994,19 +1000,22 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b merge_engine.MaskIceRings(masked_ice_rings, config_.spot_finding.ice_ring_width_Q_recipA); if (result.consensus_cell.has_value()) merge_engine.ReferenceCell(*result.consensus_cell); + // Merge-consistency filter: on floods with many spurious lattices (e.g. XFEL large cells) + // most indexed crystals do not correlate with the true structure; --min-image-cc drops the + // crystals whose per-image CC to the reference is below the limit, so the merge keeps only + // the coherent (real) ones. Set before the error model so the model, the merge and the + // reported statistics are all fitted over the same set of images. The P1 search pass keeps + // every image: its job is to find the symmetry, not to produce final intensities. + merge_engine.FilterByImageCC(!for_search + && experiment_.GetScalingSettings().GetMinCCForImage() > 0.0); merge_engine.RefineErrorModel(merge_input); if (merge_engine.ErrorModelActive()) logger.Info("Error model: a={:.3f} b={:.3f} ISa={:.1f} chi2={:.2f}", merge_engine.ErrorModelA(), merge_engine.ErrorModelB(), merge_engine.ErrorModelB() > 0 ? 1.0 / merge_engine.ErrorModelB() : 0.0, merge_engine.ErrorModelChi2()); - // Merge-consistency filter: on floods with many spurious lattices (e.g. XFEL large cells) - // most indexed crystals do not correlate with the true structure; --min-image-cc drops the - // crystals whose per-image CC to the reference is below the limit, so the merge keeps only - // the coherent (real) ones. - const bool cc_filter = !for_search && experiment_.GetScalingSettings().GetMinCCForImage() > 0.0; for (size_t i = 0; i < merge_input.size(); ++i) - merge_engine.AddImage(merge_input[i], static_cast(i), cc_filter); + merge_engine.AddImage(merge_input[i], static_cast(i)); ScaleMergeResult out; out.merged = merge_engine.ExportReflections(); diff --git a/rugnux/RugnuxCommandLine.cpp b/rugnux/RugnuxCommandLine.cpp index 751821fe..13374fa2 100644 --- a/rugnux/RugnuxCommandLine.cpp +++ b/rugnux/RugnuxCommandLine.cpp @@ -76,7 +76,10 @@ std::string RugnuxCommandLine(const ProcessConfig &config, if (azint) { const auto a = experiment.GetAzimuthalIntegrationSettings(); add("--azim-min-q", num(a.GetLowQ_recipA())); - add("--azim-max-q", num(a.GetHighQ_recipA())); + // An unset maximum Q means "to the detector edge"; emitting the resolved number would pin it + // to this run's geometry, so leave the flag out and let it resolve again. + if (const auto high_q = a.GetRequestedHighQ_recipA()) + add("--azim-max-q", num(*high_q)); add("--azim-q-spacing", num(a.GetQSpacing_recipA())); add("--azim-phi-bins", std::to_string(a.GetAzimuthalBinCount())); add("--polarization-correction", a.IsPolarizationCorrection() ? "on" : "off"); @@ -85,13 +88,18 @@ std::string RugnuxCommandLine(const ProcessConfig &config, const auto &sf = config.spot_finding; add("--spot-sigma", num(sf.signal_to_noise_threshold)); add("--spot-threshold", std::to_string(sf.photon_count_threshold)); + // The adaptive-detection default depends on the workflow (on for stills, off for rotation), so + // always emit the explicit flag rather than relying on it. + args.emplace_back(sf.adaptive_threshold ? "--adaptive-spots" : "--no-adaptive-spots"); 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)); + // Same for the spot-finding limit: absent means "as far as the detector reaches". + if (sf.high_resolution_limit.has_value()) + add("--spot-high-resolution", num(*sf.high_resolution_limit)); add("--max-spots", std::to_string(experiment.GetMaxSpotCount())); const auto idx = experiment.GetIndexingSettings(); @@ -151,8 +159,6 @@ std::string RugnuxCommandLine(const ProcessConfig &config, const auto sc = experiment.GetScalingSettings(); if (!sc.GetMergeFriedel()) args.emplace_back("-A"); - if (sc.GetRefineB()) - args.emplace_back("-B"); if (!sc.GetStillsPartialityRefine()) args.emplace_back("--simple-stills"); if (!sc.GetExpectedVarianceMerge()) diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 15b579d7..d1e1e66f 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -66,9 +66,10 @@ void print_usage() { 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. 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 << " --adaptive-spots Self-calibrating detection (DEFAULT for stills): the strong-pixel threshold comes from each image's own per-resolution-ring noise instead of the fixed --spot-threshold, so one setting adapts across datasets (no per-dataset --spot-threshold/--spot-sigma tuning). Rotation data keeps the fixed-threshold finder unless this is given." << std::endl; + std::cout << " --no-adaptive-spots Turn adaptive detection off and use the fixed --spot-threshold / --spot-sigma finder instead" << 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; + std::cout << " --spot-high-resolution High resolution limit for spot finding. If omitted, stills extend as far as the detector reaches (no resolution clipping) and rotation data keeps a 1.5 A limit." << std::endl; std::cout << " --spot-low-resolution Low resolution limit for spot finding, in A (default: 50; lower it, e.g. 24, to exclude the direct-beam halo on weakly-diffracting serial data)" << std::endl; std::cout << " --max-spots Max spot count (default: 250)" << std::endl; std::cout << " --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 dataset value)" << std::endl; @@ -99,7 +100,6 @@ void print_usage() { std::cout << " --no-scaling-corrections rot3d: disable the (default-on) decay + absorption correction surfaces fitted on the fulls after scale-fulls" << std::endl; std::cout << " --no-expected-variance-merge stills: disable the default expected-variance merge weighting (which rebuilds each weak observation's signal variance at the reflection mean to de-bias the inverse-variance merge); restores observed-sigma weighting" << std::endl; std::cout << " -A, --anomalous Anomalous mode (don't merge Friedel pairs)" << std::endl; - std::cout << " -B, --refine-bfactor Refine per image B-factor (stills only)" << std::endl; std::cout << " --scaling-high-resolution High resolution limit for scaling/merging (manual override; default: no limit)" << std::endl; std::cout << " --resolution-cutoff Automatic high-resolution cutoff for the written reflections + reported shells: cc-logistic|off (default: cc-logistic; ignored when --scaling-high-resolution is set)" << std::endl; std::cout << " --resolution-cc-target CC1/2 target defining the cc-logistic fall-off (default: 0.30)" << std::endl; @@ -124,7 +124,7 @@ void print_usage() { std::cout << " --simple-stills stills: treat every reflection as a full (p=1, single-pass scale/merge); disables the default physical partiality post-refinement" << std::endl; std::cout << " -q, --azim-q-spacing Azimuthal-integration Q bin spacing (1/A) (default: 0.01)" << std::endl; std::cout << " --azim-min-q Azimuthal-integration minimum Q (1/A)" << std::endl; - std::cout << " --azim-max-q Azimuthal-integration maximum Q (1/A)" << std::endl; + std::cout << " --azim-max-q Azimuthal-integration maximum Q (1/A). If omitted, integration extends to the highest Q the detector reaches." << std::endl; std::cout << " --azim-phi-bins Number of azimuthal (phi) bins (default: 1)" << std::endl; std::cout << " --polarization-correction Enable/disable azimuthal polarization correction" << std::endl; std::cout << " --solid-angle-correction Enable/disable azimuthal solid angle correction" << std::endl; @@ -145,6 +145,7 @@ enum { OPT_SPOT_THRESHOLD, OPT_MIN_PIX_PER_SPOT, OPT_ADAPTIVE_SPOTS, + OPT_NO_ADAPTIVE_SPOTS, OPT_SPOT_FALSE_PIXELS, OPT_SPOT_RESOLUTION, OPT_SPOT_LOW_RESOLUTION, @@ -214,7 +215,6 @@ static option long_options[] = { {"dump-observations", required_argument, nullptr, OPT_DUMP_OBSERVATIONS}, {"space-group", required_argument, nullptr, 'S'}, {"anomalous", no_argument, nullptr, 'A'}, - {"refine-bfactor", no_argument, nullptr, 'B'}, {"azint-only", no_argument, nullptr, OPT_AZINT_ONLY}, {"scale", no_argument, nullptr, OPT_SCALE}, {"no-merge", no_argument, nullptr, OPT_NO_MERGE}, @@ -253,6 +253,7 @@ static option long_options[] = { {"spot-threshold", required_argument, nullptr, OPT_SPOT_THRESHOLD}, {"min-pix-per-spot", required_argument, nullptr, OPT_MIN_PIX_PER_SPOT}, {"adaptive-spots", no_argument, nullptr, OPT_ADAPTIVE_SPOTS}, + {"no-adaptive-spots", no_argument, nullptr, OPT_NO_ADAPTIVE_SPOTS}, {"spot-false-pixels", required_argument, nullptr, OPT_SPOT_FALSE_PIXELS}, {"spot-high-resolution", required_argument, nullptr, OPT_SPOT_RESOLUTION}, {"spot-low-resolution", required_argument, nullptr, OPT_SPOT_LOW_RESOLUTION}, @@ -515,9 +516,8 @@ int main(int argc, char **argv) { float sigma_spot_finding = 3.0; int64_t photon_count_threshold_spot_finding = 10; std::optional min_pix_per_spot; // unset -> adaptive per image; a value -> fixed min-pix - bool adaptive_spots = false; + std::optional adaptive_spots; // unset -> per-workflow default (stills on, rotation off) float false_pixels_per_frame = 100.0f; - bool refine_bfactor = false; std::string ref_mtz; std::string ref_column; std::string model_pdb; // --model: PDB to validate merged intensities against (R-free + maps) @@ -537,7 +537,7 @@ int main(int argc, char **argv) { IndexingAlgorithmEnum indexing_algorithm = IndexingAlgorithmEnum::Auto; GeomRefinementAlgorithmEnum refinement_algorithm = GeomRefinementAlgorithmEnum::BeamCenter; - float d_min_spot_finding = 1.5; + std::optional d_min_spot_finding; // unset -> as far as the detector reaches float d_max_spot_finding = 0; // 0 = keep the SpotFindingSettings default (50 A) std::optional d_min_scale_merge; std::optional resolution_cutoff_method; // --resolution-cutoff cc-logistic|off @@ -556,7 +556,7 @@ int main(int argc, char **argv) { int opt; int option_index = 0; - const char *short_opts = "vo:N:s:e:t:R::X:C:z:FABS:r:q:"; + const char *short_opts = "vo:N:s:e:t:R::X:C:z:FAS:r:q:"; while ((opt = getopt_long(argc, argv, short_opts, long_options, &option_index)) != -1) { switch (opt) { @@ -715,9 +715,6 @@ int main(int argc, char **argv) { case 'A': anomalous_mode = true; break; - case 'B': - refine_bfactor = true; - break; case 'S': { // Accept a space-group number ("92") or a Hermann-Mauguin symbol ("P43212", "P 43 21 2"). char *end = nullptr; @@ -750,6 +747,10 @@ int main(int argc, char **argv) { adaptive_spots = true; logger.Info("Adaptive (self-calibrating) spot detection enabled"); break; + case OPT_NO_ADAPTIVE_SPOTS: + adaptive_spots = false; + logger.Info("Adaptive spot detection off: using the fixed --spot-threshold / --spot-sigma finder"); + break; case OPT_SPOT_FALSE_PIXELS: false_pixels_per_frame = parse_number_arg(optarg, "--spot-false-pixels", logger, 1.0f); adaptive_spots = true; @@ -759,10 +760,20 @@ int main(int argc, char **argv) { d_max_spot_finding = parse_number_arg(optarg, "--spot-low-resolution", logger, 0.0f); logger.Info("Low resolution limit for spot finding set to {:.1f} A", d_max_spot_finding); break; - case OPT_SPOT_RESOLUTION: - d_min_spot_finding = parse_number_arg(optarg, "--spot-high-resolution", logger, 0.0f); - logger.Info("High resolution limit for spot finding set to {:.2f} A", d_min_spot_finding); + case OPT_SPOT_RESOLUTION: { + // 0 has always meant "no limit" for this setting; keep that, but express it as the unset + // optional the rest of the code understands. Passing the 0 through instead reached + // ResolutionShells (via the spot plot), which rejects a zero d_min and threw away every image. + const auto d_min = parse_number_arg(optarg, "--spot-high-resolution", logger, 0.0f); + if (d_min > 0.0f) { + d_min_spot_finding = d_min; + logger.Info("High resolution limit for spot finding set to {:.2f} A", d_min); + } else { + d_min_spot_finding.reset(); + logger.Info("No high resolution limit for spot finding: as far as the detector reaches"); + } break; + } case OPT_MAX_SPOTS: max_spot_count_override = parse_number_arg(optarg, "--max-spots", logger, 1); break; @@ -1074,7 +1085,6 @@ int main(int argc, char **argv) { if (resolution_cc_target) scaling_settings.ResolutionCCTarget(*resolution_cc_target); if (report_shell_count) scaling_settings.ReportShellCount(*report_shell_count); scaling_settings.MergeFriedel(!anomalous_mode); - scaling_settings.RefineB(refine_bfactor); scaling_settings.MinPartiality(min_partiality); scaling_settings.MinCapturedFraction(min_captured_fraction_arg.value_or( (experiment.GetGoniometer().has_value() && !force_still) ? 0.7 : 0.0)); @@ -1112,16 +1122,16 @@ int main(int argc, char **argv) { // Rotation (rot3d): the dedicated RotationScaleMerge does the whole self-scale -> 3D combine -> // merge, including the default-on decay + absorption correction surfaces. It does not support - // external-reference scaling, the stills -B (per-image B-factor) or wedge refinement. + // external-reference scaling or wedge refinement. // Everything else (stills, reference scaling) uses ScaleOnTheFly + MergeOnTheFly. const bool is_rotation = experiment.IsRotationIndexing(); if (is_rotation) { - if (!reference_data.empty() || experiment.GetScalingSettings().GetRefineB() + if (!reference_data.empty() || experiment.GetRefineRotationWedgeInScaling() || experiment.GetScalingSettings().GetRotationWedgeForScaling().has_value()) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Rotation scaling/merging (RotationScaleMerge) does not support reference " - "scaling, B-factor refinement or wedge refinement"); + "scaling or wedge refinement"); RotationScaleMerge rsm(experiment, reflections, experiment.GetUnitCell(), scaling_iter, 0.0f, nthreads, logger); rsm.Ingest(); @@ -1146,6 +1156,8 @@ int main(int argc, char **argv) { } MergeOnTheFly merge_engine(experiment); merge_engine.ReferenceCell(experiment.GetUnitCell()); + // --min-image-cc has to hold for the merge itself, not only for the reported statistics. + merge_engine.FilterByImageCC(experiment.GetScalingSettings().GetMinCCForImage() > 0.0); // Fit the (a, b) error model from symmetry-mate scatter before merging, exactly as the full // pipeline does (Rugnux.cpp). Without this the offline --scale merge would use the identity // model and produce much worse stills intensities (no (b*I)^2 systematic term, no sigma floor). @@ -1262,7 +1274,7 @@ int main(int argc, char **argv) { AzimuthalIntegrationSettings azint_settings = experiment.GetAzimuthalIntegrationSettings(); if (min_q || max_q) azint_settings.QRange_recipA(min_q.value_or(azint_settings.GetLowQ_recipA()), - max_q.value_or(azint_settings.GetHighQ_recipA())); + max_q ? max_q : azint_settings.GetRequestedHighQ_recipA()); if (q_spacing) azint_settings.QSpacing_recipA(q_spacing.value()); if (azimuthal_bins) @@ -1398,7 +1410,6 @@ int main(int argc, char **argv) { if (resolution_cc_target) scaling_settings.ResolutionCCTarget(*resolution_cc_target); if (report_shell_count) scaling_settings.ReportShellCount(*report_shell_count); scaling_settings.MergeFriedel(!anomalous_mode); - scaling_settings.RefineB(refine_bfactor); scaling_settings.MinPartiality(min_partiality); // Drop edge-of-sweep truncated fulls (rocking curve captured < this fraction) from the rot3d combine. // Defaults ON (0.7) for rotation - removes the low-capture fulls that inflate low-res R-meas and @@ -1468,14 +1479,25 @@ int main(int argc, char **argv) { SpotFindingSettings spot_settings; spot_settings.enable = true; spot_settings.indexing = true; - spot_settings.high_resolution_limit = d_min_spot_finding; spot_settings.signal_to_noise_threshold = sigma_spot_finding; spot_settings.photon_count_threshold = photon_count_threshold_spot_finding; + // Detection defaults differ by workflow; each is overridden by its flag, which always wins. + // - min-pix: choosing it per image (unset) only means something where each frame is indexed on its + // own. Rotation indexing builds ONE lattice from all frames, so it keeps the fixed value. + // - adaptive detection: a clear win across the stills battery, but on the 33-crystal rotation + // battery it deterministically breaks three (a lost space group, a halved indexing rate, a + // collapsed merge) while helping four, so rotation keeps the fixed-threshold finder. + // - high-resolution limit: unset means "as far as the detector reaches", which is what the stills + // path wants; on rotation the extra high-resolution spots cost indexing (measured 100.0 -> 96.8% + // on a strong rotation set), so rotation keeps the historical limit. spot_settings.min_pix_per_spot = min_pix_per_spot; - spot_settings.adaptive_threshold = adaptive_spots; + if (rotation_indexing && !spot_settings.min_pix_per_spot.has_value()) + spot_settings.min_pix_per_spot = 2; + if (rotation_indexing && !d_min_spot_finding.has_value()) + d_min_spot_finding = 1.5f; + spot_settings.adaptive_threshold = adaptive_spots.value_or(!rotation_indexing); + spot_settings.high_resolution_limit = d_min_spot_finding; spot_settings.false_pixels_per_frame = false_pixels_per_frame; - if (d_min_spot_finding > 0.0f) - spot_settings.high_resolution_limit = d_min_spot_finding; if (d_max_spot_finding > 0.0f) spot_settings.low_resolution_limit = d_max_spot_finding; diff --git a/tests/AdaptiveThresholdTest.cpp b/tests/AdaptiveThresholdTest.cpp new file mode 100644 index 00000000..ecd13da6 --- /dev/null +++ b/tests/AdaptiveThresholdTest.cpp @@ -0,0 +1,100 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#include + +#include "../image_analysis/spot_finding/AdaptiveThreshold.h" + +using namespace adaptive_threshold; + +namespace { + // Poisson upper tail P(X >= k) for mean mu, summed directly - an independent reference for the + // threshold's defining property. + double PoissonUpperTail(double mu, int k) { + if (k <= 0) + return 1.0; + double pmf = std::exp(-mu); + double cdf = pmf; + for (int i = 1; i < k; i++) { + pmf *= mu / i; + cdf += pmf; + } + return std::max(0.0, 1.0 - cdf); + } +} + +TEST_CASE("AdaptiveThreshold_NormalQuantile", "[SpotFinding]") { + // Textbook values of the inverse standard-normal CDF. + CHECK(NormalQuantile(0.5) == Catch::Approx(0.0).margin(1e-9)); + CHECK(NormalQuantile(0.975) == Catch::Approx(1.959964).margin(1e-5)); + CHECK(NormalQuantile(0.99) == Catch::Approx(2.326348).margin(1e-5)); + CHECK(NormalQuantile(1.0 - 1e-6) == Catch::Approx(4.753424).margin(1e-4)); + + // Symmetric about 0.5, and monotonically increasing. + for (const double p: {1e-8, 1e-4, 0.01, 0.2, 0.45}) + CHECK(NormalQuantile(1.0 - p) == Catch::Approx(-NormalQuantile(p)).margin(1e-6)); + CHECK(NormalQuantile(0.6) > NormalQuantile(0.55)); + CHECK(NormalQuantile(1e-3) < NormalQuantile(1e-2)); + + // Degenerate arguments stay finite: the finders divide a tolerated-false-pixel count by the pixel + // count, so p can legitimately arrive at the very edge of (0, 1). + CHECK(std::isfinite(NormalQuantile(0.0))); + CHECK(std::isfinite(NormalQuantile(1.0))); + CHECK(NormalQuantile(0.0) < 0.0); + CHECK(NormalQuantile(1.0) > 0.0); +} + +TEST_CASE("AdaptiveThreshold_PoissonThreshold", "[SpotFinding]") { + const double p = 1e-5; + 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}) { + const int thr = static_cast(PoissonThreshold(mu, p, z)); + CHECK(PoissonUpperTail(mu, thr) <= p); + CHECK(PoissonUpperTail(mu, thr - 1) > p); + } + + // Non-decreasing in the background level. + float prev = 0.0f; + for (const double mu: {1e-6, 0.01, 0.1, 0.5, 1.0, 2.0, 5.0, 20.0, 45.0}) { + const float thr = PoissonThreshold(mu, p, z); + CHECK(thr >= prev); + 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)); + + // A tighter operating point (smaller p) can only raise the threshold. + CHECK(PoissonThreshold(5.0, 1e-8, static_cast(NormalQuantile(1.0 - 1e-8))) + >= PoissonThreshold(5.0, 1e-2, static_cast(NormalQuantile(1.0 - 1e-2)))); +} + +TEST_CASE("AdaptiveThreshold_RingThreshold", "[SpotFinding]") { + const double p = 1e-5; + const float z = static_cast(NormalQuantile(1.0 - p)); + + // Never below the read-noise-aware Gaussian arm, which is what keeps an empty ring's threshold + // off zero - a per-ring sigma alone would collapse there and flood the frame with noise spots. + for (const float mean: {0.0f, 0.5f, 5.0f, 50.0f}) { + for (const float sigma: {0.0f, 1.0f, 7.0f}) { + const float gauss = mean + z * std::sqrt(sigma * sigma + READ * READ); + CHECK(RingThreshold(mean, sigma, p, z) >= Catch::Approx(gauss).epsilon(1e-6)); + } + } + CHECK(RingThreshold(0.0f, 0.0f, p, z) >= z * READ); + + // Non-decreasing in the background mean and in the background scatter. + CHECK(RingThreshold(20.0f, 4.0f, p, z) > RingThreshold(2.0f, 4.0f, p, z)); + CHECK(RingThreshold(5.0f, 9.0f, p, z) > RingThreshold(5.0f, 1.0f, p, z)); + + // Where the background is countable and quiet, Poisson significance is the binding arm: a ring + // with mean 1 and no measured scatter must still demand several photons. + CHECK(RingThreshold(1.0f, 0.0f, p, z) > 1.0f + z * READ); + + // A ring whose scatter is far above Poisson (flat-field / read excess) is set by the Gaussian arm. + CHECK(RingThreshold(10.0f, 30.0f, p, z) == Catch::Approx(10.0f + z * std::sqrt(900.0f + READ * READ)).epsilon(1e-6)); +} diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 0dade2aa..9145a49b 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -6,6 +6,7 @@ ADD_EXECUTABLE(jfjoch_test CheckImageOutput.h FPGAIntegrationTest.cpp StrongPixelSetTest.cpp + AdaptiveThresholdTest.cpp ZSTDCompressorTest.cpp FrameTransformationTest.cpp HDF5WritingTest.cpp PedestalCalcTest.cpp @@ -64,6 +65,7 @@ ADD_EXECUTABLE(jfjoch_test ImageSpotFinderCPUTest.cpp ImageSpotFinderGPUTest.cpp AdaptiveSpotFinderGPUTest.cpp + AdaptiveThresholdTest.cpp CalcBraggPredictionTest.cpp SpotUtilsTest.cpp LatticeSearchTest.cpp diff --git a/tests/JFJochReaderTest.cpp b/tests/JFJochReaderTest.cpp index 68243fcf..e84d6a25 100644 --- a/tests/JFJochReaderTest.cpp +++ b/tests/JFJochReaderTest.cpp @@ -985,6 +985,9 @@ TEST_CASE("JFJochReader_Azint", "[HDF5][Full]") { AzimuthalIntegrationSettings azint_settings; azint_settings.AzimuthalBinCount(4); x.ImportAzimuthalIntegrationSettings(azint_settings); + // The high-q limit is unset, i.e. "as far as the detector reaches", so read the settings back from + // the experiment, where that has been resolved against the geometry - that is what the bins are. + azint_settings = x.GetAzimuthalIntegrationSettings(); std::vector image(x.GetPixelsNum()); @@ -1544,6 +1547,9 @@ TEST_CASE("JFJochReader_NXmxIntegrated", "[HDF5][Full]") { AzimuthalIntegrationSettings azint_settings; azint_settings.AzimuthalBinCount(4); x.ImportAzimuthalIntegrationSettings(azint_settings); + // The high-q limit is unset, i.e. "as far as the detector reaches", so read the settings back from + // the experiment, where that has been resolved against the geometry - that is what the bins are. + azint_settings = x.GetAzimuthalIntegrationSettings(); std::vector image(x.GetPixelsNum(), 0); image[0] = UINT16_MAX; diff --git a/tests/XDSPluginTest.cpp b/tests/XDSPluginTest.cpp index 926893bf..bd23b444 100644 --- a/tests/XDSPluginTest.cpp +++ b/tests/XDSPluginTest.cpp @@ -191,6 +191,8 @@ TEST_CASE("XDSPlugin_GetData_Integrated", "[HDF5][XDS][Plugin]") { AzimuthalIntegrationSettings azint_settings; azint_settings.AzimuthalBinCount(4); x.ImportAzimuthalIntegrationSettings(azint_settings); + // Unset high q = "as far as the detector reaches"; the experiment resolves it against the geometry. + azint_settings = x.GetAzimuthalIntegrationSettings(); std::vector image(x.GetPixelsNum(), 0); image[0] = UINT16_MAX; diff --git a/update_version.sh b/update_version.sh index e29f57d0..903c529b 100644 --- a/update_version.sh +++ b/update_version.sh @@ -4,6 +4,11 @@ # Copyright (2019-2024) Paul Scherrer Institute # +# Stop at the first failing command (this script rewrites generated code and version +# strings in place, so a silently skipped step leaves the tree half-updated). +set -e +set -o pipefail + VERSION=$( docs/THIRD_PARTY_NOTICES.md -git add broker/gen/model/*.cpp broker/gen/model/*.h frontend/src/openapi/models/*.ts docs/python_client/*.md docs/python_client/docs/*.md docs/THIRD_PARTY_NOTICES.md +git add broker/gen/model/*.cpp broker/gen/model/*.h frontend/src/client docs/python_client/*.md docs/python_client/docs/*.md docs/THIRD_PARTY_NOTICES.md sed -i "s,release =.*,release = \'$VERSION\'," docs/conf.py -sed -i "s,\"version\":.*,\"version\": \"$VERSION\"\,," frontend/package.json cd fpga diff --git a/viewer/JFJochHttpReader.cpp b/viewer/JFJochHttpReader.cpp index 4385297d..95210bf2 100644 --- a/viewer/JFJochHttpReader.cpp +++ b/viewer/JFJochHttpReader.cpp @@ -281,7 +281,6 @@ std::shared_ptr JFJochHttpReader::UpdateDataset_i() { dataset->integrated_reflections = GetPlot_i("integrated_reflections"); dataset->image_scale_factor = GetPlot_i("image_scale_factor"); dataset->image_scale_cc = GetPlot_i("image_scale_cc"); - dataset->image_scale_b = GetPlot_i("image_scale_b"); if (msg->start_message->goniometer) dataset->experiment.Goniometer(msg->start_message->goniometer); diff --git a/viewer/JFJochProcessController.cpp b/viewer/JFJochProcessController.cpp index 8b7c1016..8301737a 100644 --- a/viewer/JFJochProcessController.cpp +++ b/viewer/JFJochProcessController.cpp @@ -119,7 +119,6 @@ void JFJochProcessController::OnImageProcessed(const DataMessage &msg) { if (msg.integrated_reflections) put(d.integrated_reflections, *msg.integrated_reflections); if (msg.image_scale_factor) put(d.image_scale_factor, *msg.image_scale_factor); if (msg.image_scale_cc) put(d.image_scale_cc, *msg.image_scale_cc); - if (msg.image_scale_b_factor) put(d.image_scale_b, *msg.image_scale_b_factor); // Throttle to ~4 Hz so the GUI plots refresh smoothly without flooding the event queue. const auto now = std::chrono::steady_clock::now(); diff --git a/viewer/JFJochViewerWindow.cpp b/viewer/JFJochViewerWindow.cpp index 56ea0a8f..3501c328 100644 --- a/viewer/JFJochViewerWindow.cpp +++ b/viewer/JFJochViewerWindow.cpp @@ -85,7 +85,12 @@ JFJochViewerWindow::JFJochViewerWindow(QWidget *parent, bool dbus, const QString SpotFindingSettings spot_finding_settings = DiffractionExperiment::DefaultDataProcessingSettings(); - spot_finding_settings.high_resolution_limit = 1.5; + // The interactive viewer is an offline tool looking at one dataset at a time, so it starts from the + // self-calibrating settings rather than numbers the user would have to tune: the per-ring adaptive + // threshold, min-pix chosen per image (std::nullopt), and detection out to the detector edge + // (high_resolution_limit unset). All three are switchable in the settings dock. + spot_finding_settings.adaptive_threshold = true; + spot_finding_settings.min_pix_per_spot = std::nullopt; spot_finding_settings.indexing = true; IndexingSettings indexing_settings; diff --git a/viewer/image_viewer/JFJochDiffractionImage.cpp b/viewer/image_viewer/JFJochDiffractionImage.cpp index 267415c5..44c6085f 100644 --- a/viewer/image_viewer/JFJochDiffractionImage.cpp +++ b/viewer/image_viewer/JFJochDiffractionImage.cpp @@ -1030,11 +1030,6 @@ void JFJochDiffractionImage::DrawResolutionText() { viewport()->update(dirty); } -void JFJochDiffractionImage::beforeOverlayCleared() { - // The resolution readout is painted in drawForeground(), not held as a scene item, so - // clearing the overlay (or the whole scene) cannot leave a dangling pointer behind. -} - void JFJochDiffractionImage::leaveEvent(QEvent *event) { // Mouse left the view: clear hover resolution and hide text if (std::isfinite(hover_resolution)) { diff --git a/viewer/image_viewer/JFJochDiffractionImage.h b/viewer/image_viewer/JFJochDiffractionImage.h index f3998be1..f4cc0ad8 100644 --- a/viewer/image_viewer/JFJochDiffractionImage.h +++ b/viewer/image_viewer/JFJochDiffractionImage.h @@ -87,7 +87,6 @@ private: void DrawCross(float x, float y, float size, float width, float z = 1); void UpdateForeground(); - void beforeOverlayCleared() override; void leaveEvent(QEvent *event) override; std::shared_ptr image; diff --git a/viewer/image_viewer/JFJochFollowerImage.cpp b/viewer/image_viewer/JFJochFollowerImage.cpp index 95761db5..fd113cb9 100644 --- a/viewer/image_viewer/JFJochFollowerImage.cpp +++ b/viewer/image_viewer/JFJochFollowerImage.cpp @@ -33,7 +33,6 @@ JFJochFollowerImage::JFJochFollowerImage(QWidget *parent) : QGraphicsView(parent void JFJochFollowerImage::SetFrame(std::shared_ptr frame) { const bool same_pointer = (frame_ == frame); - const bool same_size = frame_ && frame && frame_->size() == frame->size(); frame_ = std::move(frame); if (!frame_ || frame_->isNull()) { @@ -49,9 +48,12 @@ void JFJochFollowerImage::SetFrame(std::shared_ptr frame) { scene()->addItem(item_); } - if (!same_size) { + // The producer reassigns its buffer in place when the frame size changes, so the pointer says + // nothing about the dimensions - compare against the size we last saw. + if (frame_->size() != frame_size_) { + frame_size_ = frame_->size(); item_->refresh(); - scene()->setSceneRect(0, 0, frame_->width(), frame_->height()); + scene()->setSceneRect(0, 0, frame_size_.width(), frame_size_.height()); } viewport()->update(); diff --git a/viewer/image_viewer/JFJochFollowerImage.h b/viewer/image_viewer/JFJochFollowerImage.h index 90bd5f5d..bcafb990 100644 --- a/viewer/image_viewer/JFJochFollowerImage.h +++ b/viewer/image_viewer/JFJochFollowerImage.h @@ -30,6 +30,7 @@ class JFJochFollowerImage : public QGraphicsView { JFJochImageItem *item_ = nullptr; std::shared_ptr frame_; + QSize frame_size_; // size the scene rect was set from std::shared_ptr values_; double zoom_ = 12.0; diff --git a/viewer/image_viewer/JFJochImage.cpp b/viewer/image_viewer/JFJochImage.cpp index 9fb40eba..25077bae 100644 --- a/viewer/image_viewer/JFJochImage.cpp +++ b/viewer/image_viewer/JFJochImage.cpp @@ -94,6 +94,12 @@ void JFJochImage::ScheduleHoverUpdate(const QPointF &scenePos, Qt::KeyboardModif hover_tail_timer_->start(kHoverIntervalMs); } +void JFJochImage::leaveEvent(QEvent *event) { + // A tail update that fires now would report a position the pointer has already left + hover_tail_timer_->stop(); + QGraphicsView::leaveEvent(event); +} + void JFJochImage::ScheduleRenderImage() { if (render_pending_) return; @@ -900,8 +906,6 @@ void JFJochImage::resetScenePointers() { void JFJochImage::updateOverlay() { if (!scene() || W * H <= 0) return; - beforeOverlayCleared(); - // Remove only overlay items, keep the image item persistent for (auto *item : overlay_items_) scene()->removeItem(item); @@ -975,8 +979,6 @@ void JFJochImage::adjustForeground(bool input) { m_adjustForegroundWithWheel = input; } -void JFJochImage::beforeOverlayCleared() {} - double JFJochImage::GetScaleFactor() const { return scale_factor; } diff --git a/viewer/image_viewer/JFJochImage.h b/viewer/image_viewer/JFJochImage.h index 23f8ddba..b3ead6b1 100644 --- a/viewer/image_viewer/JFJochImage.h +++ b/viewer/image_viewer/JFJochImage.h @@ -95,7 +95,6 @@ protected: // Only the view that owns detector counts offers a region of interest; for the others a // shift-drag would draw a box that means nothing. [[nodiscard]] virtual bool AllowROI() const { return false; } - virtual void beforeOverlayCleared(); bool show_saturation = false; @@ -178,6 +177,7 @@ protected: static constexpr int kHoverIntervalMs = 66; void ScheduleHoverUpdate(const QPointF &scenePos, Qt::KeyboardModifiers modifiers); void UpdateHover(); + void leaveEvent(QEvent *event) override; QPointF hover_scene_pos_; Qt::KeyboardModifiers hover_modifiers_ = Qt::NoModifier; QElapsedTimer hover_rate_; diff --git a/viewer/widgets/JFJochViewerSettingsDock.cpp b/viewer/widgets/JFJochViewerSettingsDock.cpp index 3203b55b..25c3c385 100644 --- a/viewer/widgets/JFJochViewerSettingsDock.cpp +++ b/viewer/widgets/JFJochViewerSettingsDock.cpp @@ -55,7 +55,9 @@ JFJochViewerSettingsDock::JFJochViewerSettingsDock(const SpotFindingSettings &sp const BraggIntegrationSettings &bragg, const ScalingSettings &scaling, QWidget *parent) - : QWidget(parent), spot_(spot), indexing_(indexing), azint_(azint), bragg_(bragg), scaling_(scaling) { + : QWidget(parent), spot_(spot), indexing_(indexing), azint_(azint), bragg_(bragg), scaling_(scaling), + adaptive_min_pix_(!spot.min_pix_per_spot.has_value()), + min_pix_value_(spot.min_pix_per_spot.value_or(2)) { auto *layout = new QVBoxLayout(this); @@ -228,19 +230,28 @@ QWidget *JFJochViewerSettingsDock::BuildMXPage() { adaptive->setChecked(spot_.adaptive_threshold); adaptive->setToolTip("Self-calibrating detection: the strong-pixel threshold is derived from each " "image's own per-resolution-ring noise. The signal/noise and photon-count " - "settings are not used while this is on."); + "settings are not used while this is on. Best on stills; on rotation data the " + "fixed-threshold finder is the safer choice, which is what rugnux defaults to " + "there."); 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.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."); + highResSpot->setValue(spot_.high_resolution_limit.value_or(1.5f)); + auto *autoHighResSpot = new QCheckBox("To detector edge", page); + autoHighResSpot->setChecked(!spot_.high_resolution_limit.has_value()); + autoHighResSpot->setToolTip("Find spots as far as the detector reaches, instead of clipping the " + "detection at a fixed resolution. The high-resolution value is not used " + "while this is on."); + auto *minPix = new NumberLineEdit(1.0f, 50.0f, static_cast(min_pix_value_), 0, "px", page); + auto *adaptiveMinPix = new QCheckBox("Adaptive min-pix (stills)", page); + adaptiveMinPix->setChecked(adaptive_min_pix_); + adaptiveMinPix->setToolTip("Choose the minimum pixels/spot per image: index the frame at min-pix " + "3/2/1 and keep whichever maximises indexed count x indexed fraction. " + "Stills only - a rotation dataset builds one lattice from all frames and " + "always uses the fixed min-pixels/spot value below."); 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("", autoHighResSpot); spot->addRow("High resolution [Å]", highResSpot); spot->addRow("", adaptiveMinPix); spot->addRow("Min pixels/spot", minPix); @@ -255,7 +266,7 @@ QWidget *JFJochViewerSettingsDock::BuildMXPage() { connect(highResSpot, &SliderPlusBox::valueChanged, this, [this](double v) { spot_.high_resolution_limit = static_cast(v); EmitSpotFinding(); }); connect(minPix, &NumberLineEdit::newValue, this, [this, minPix] { - spot_.min_pix_per_spot = std::llround(minPix->value()); EmitSpotFinding(); }); + min_pix_value_ = std::llround(minPix->value()); EmitSpotFinding(); }); connect(maxSpots, &NumberLineEdit::newValue, this, [this, maxSpots] { max_spots_ = std::llround(maxSpots->value()); EmitSpotFinding(); }); // The adaptive finder sets its own threshold from each image's noise, so the signal/noise and @@ -270,16 +281,23 @@ 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) { + // An unset high-resolution limit (std::nullopt) means "to the detector edge", so the value is unused + // while that is on - grey it out, as for the other automatic settings. + auto syncHighResEnabled = [highResSpot](bool auto_on) { highResSpot->setEnabled(!auto_on); }; + syncHighResEnabled(!spot_.high_resolution_limit.has_value()); + connect(autoHighResSpot, &QCheckBox::toggled, this, [this, highResSpot, syncHighResEnabled](bool on) { if (on) - spot_.min_pix_per_spot = std::nullopt; + spot_.high_resolution_limit = std::nullopt; else - spot_.min_pix_per_spot = std::llround(minPix->value()); - syncMinPixEnabled(on); + spot_.high_resolution_limit = static_cast(highResSpot->value()); + syncHighResEnabled(on); + EmitSpotFinding(); + }); + + // The fixed min-pixels/spot field stays live even with adaptive min-pix on: it is what a rotation + // dataset uses, and what stills fall back to when adaptive is switched off. + connect(adaptiveMinPix, &QCheckBox::toggled, this, [this](bool on) { + adaptive_min_pix_ = on; EmitSpotFinding(); }); @@ -395,6 +413,10 @@ QWidget *JFJochViewerSettingsDock::BuildAzIntPage() { lowQ->setValue(azint_.GetLowQ_recipA()); auto *highQ = new SliderPlusBox(2e-5, 10.0, 0.001, 4, page); highQ->setValue(azint_.GetHighQ_recipA()); + auto *autoHighQ = new QCheckBox("To detector edge", page); + autoHighQ->setChecked(!azint_.GetRequestedHighQ_recipA().has_value()); + autoHighQ->setToolTip("Integrate out to the highest Q the detector reaches. The high-Q value is not " + "used while this is on."); auto *spacing = new SliderPlusBox(1e-5, 1.0, 0.001, 5, page, SliderPlusBox::ScaleType::Logarithmic); spacing->setValue(azint_.GetQSpacing_recipA()); auto *azimBins = new QComboBox(page); @@ -402,6 +424,7 @@ QWidget *JFJochViewerSettingsDock::BuildAzIntPage() { azimBins->addItem(QString::number(b), b); azimBins->setCurrentIndex(azimBins->findData(azint_.GetAzimuthalBinCount())); az->addRow("Low Q [Å⁻¹]", lowQ); + az->addRow("", autoHighQ); az->addRow("High Q [Å⁻¹]", highQ); az->addRow("Q spacing [Å⁻¹]", spacing); az->addRow("Azimuthal bins", azimBins); @@ -410,7 +433,9 @@ QWidget *JFJochViewerSettingsDock::BuildAzIntPage() { layout->addWidget(azSection); auto emitAz = [=, this] { - azint_.QRange_recipA(static_cast(lowQ->value()), static_cast(highQ->value())); + azint_.QRange_recipA(static_cast(lowQ->value()), + autoHighQ->isChecked() ? std::nullopt + : std::optional(highQ->value())); azint_.QSpacing_recipA(static_cast(spacing->value())); azint_.AzimuthalBinCount(azimBins->currentData().toInt()); emit azintChanged(azint_); @@ -419,6 +444,11 @@ QWidget *JFJochViewerSettingsDock::BuildAzIntPage() { connect(highQ, &SliderPlusBox::valueChanged, this, [emitAz] { emitAz(); }); connect(spacing, &SliderPlusBox::valueChanged, this, [emitAz] { emitAz(); }); connect(azimBins, &QComboBox::currentIndexChanged, this, [emitAz] { emitAz(); }); + connect(autoHighQ, &QCheckBox::toggled, this, [emitAz, highQ](bool on) { + highQ->setEnabled(!on); + emitAz(); + }); + highQ->setEnabled(!autoHighQ->isChecked()); // Powder calibration (calibrant rings + geometry refinement) - reuse the existing widget. auto *powderSection = new CollapsibleSection("Powder calibration", page); @@ -436,7 +466,16 @@ QWidget *JFJochViewerSettingsDock::BuildAzIntPage() { return page; } +void JFJochViewerSettingsDock::SyncMinPix() { + // Per-image min-pix indexes each frame on its own, which only means something for stills; rotation + // indexing builds one lattice from all frames, so it keeps the fixed value. + spot_.min_pix_per_spot = (adaptive_min_pix_ && !indexing_.GetRotationIndexing()) + ? std::optional() + : std::optional(min_pix_value_); +} + void JFJochViewerSettingsDock::EmitSpotFinding() { + SyncMinPix(); emit spotFindingChanged(spot_, indexing_, max_spots_); } @@ -528,8 +567,6 @@ QWidget *JFJochViewerSettingsDock::BuildScalingSection() { auto *friedel = new QCheckBox("Merge Friedel pairs", this); friedel->setChecked(scaling_.GetMergeFriedel()); - auto *refineB = new QCheckBox("Refine B-factor", this); - refineB->setChecked(scaling_.GetRefineB()); auto *corrections = new QCheckBox("Correction surfaces (decay + absorption)", this); corrections->setChecked(scaling_.GetCorrectionSurfaces()); corrections->setToolTip("Rotation only: fit a radiation-damage decay and a goniometer-frame absorption " @@ -547,7 +584,6 @@ QWidget *JFJochViewerSettingsDock::BuildScalingSection() { highRes->setEnabled(limitRes->isChecked()); form->addRow("", friedel); - form->addRow("", refineB); form->addRow("", corrections); form->addRow("", partRefine); // Compact, and aligned with the checkboxes above: the limit checkbox + value sit together in the @@ -561,7 +597,6 @@ QWidget *JFJochViewerSettingsDock::BuildScalingSection() { auto emitScaling = [=, this] { scaling_.MergeFriedel(friedel->isChecked()); - scaling_.RefineB(refineB->isChecked()); scaling_.CorrectionSurfaces(corrections->isChecked()); scaling_.StillsPartialityRefine(partRefine->isChecked()); scaling_.HighResolutionLimit_A(limitRes->isChecked() @@ -569,7 +604,6 @@ QWidget *JFJochViewerSettingsDock::BuildScalingSection() { emit scalingChanged(scaling_); }; connect(friedel, &QCheckBox::toggled, this, [emitScaling] { emitScaling(); }); - connect(refineB, &QCheckBox::toggled, this, [emitScaling] { emitScaling(); }); connect(corrections, &QCheckBox::toggled, this, [emitScaling] { emitScaling(); }); connect(partRefine, &QCheckBox::toggled, this, [emitScaling] { emitScaling(); }); connect(limitRes, &QCheckBox::toggled, this, [emitScaling, highRes](bool on) { diff --git a/viewer/widgets/JFJochViewerSettingsDock.h b/viewer/widgets/JFJochViewerSettingsDock.h index d4e36b2a..c95b67d1 100644 --- a/viewer/widgets/JFJochViewerSettingsDock.h +++ b/viewer/widgets/JFJochViewerSettingsDock.h @@ -65,6 +65,11 @@ private: DiffractionExperiment experiment_; bool have_experiment_ = false; // geometry edits only take effect once a dataset is loaded int64_t max_spots_ = 1000; + // Min-pix is kept as the pair the UI shows - the fixed value AND whether to choose it per image - + // because spot_.min_pix_per_spot can only hold one of the two (unset = per image). SyncMinPix() + // turns the pair into that field, and per-image only ever applies to stills. + bool adaptive_min_pix_ = false; + int64_t min_pix_value_ = 2; bool azint_mode_ = false; // false = MX page, true = AzInt page (drives "Analyze dataset") // "Analyze dataset" hero button, disabled while a live HTTP source is connected. @@ -107,6 +112,7 @@ private: QWidget *BuildScalingSection(); QWidget *BuildReferenceSection(); QWidget *BuildAzIntPage(); + void SyncMinPix(); void EmitSpotFinding(); void EmitExperiment(); void ApplyProcessingMode(); // "Process as stills" -> indexing + scaling rotation/stills mode diff --git a/viewer/windows/JFJochViewerImageListWindow.cpp b/viewer/windows/JFJochViewerImageListWindow.cpp index abd2dc5c..c46ee440 100644 --- a/viewer/windows/JFJochViewerImageListWindow.cpp +++ b/viewer/windows/JFJochViewerImageListWindow.cpp @@ -51,7 +51,7 @@ JFJochViewerImageListWindow::JFJochViewerImageListWindow(QWidget *parent) : JFJo void JFJochViewerImageListWindow::setupTableModel() { - tableModel->setColumnCount(9); + tableModel->setColumnCount(8); tableModel->setHeaderData(0, Qt::Horizontal, "#"); tableModel->setHeaderData(1, Qt::Horizontal, "Bkg"); tableModel->setHeaderData(2, Qt::Horizontal, "Index"); @@ -60,7 +60,6 @@ void JFJochViewerImageListWindow::setupTableModel() tableModel->setHeaderData(5, Qt::Horizontal, "Max val"); tableModel->setHeaderData(6, Qt::Horizontal, "Scale factor"); tableModel->setHeaderData(7, Qt::Horizontal, "Scale CC"); -tableModel->setHeaderData(8, Qt::Horizontal, "Scale B [A^2]"); } void JFJochViewerImageListWindow::addDataRow(int imageNumber, double backgroundEstimate, @@ -69,8 +68,7 @@ void JFJochViewerImageListWindow::addDataRow(int imageNumber, double backgroundE double resolutionEstimate, int64_t max_value, double image_scale_factor, - double image_scale_cc_percent, - double image_scale_b) { + double image_scale_cc_percent) { QList rowItems; QStandardItem *imageItem = new QStandardItem(); @@ -128,16 +126,6 @@ void JFJochViewerImageListWindow::addDataRow(int imageNumber, double backgroundE } rowItems.append(scaleCCItem); - QStandardItem *scaleBItem = new QStandardItem(); - if (std::isfinite(image_scale_b)) { - scaleBItem->setData(QString::number(image_scale_b, 'f', 2), Qt::DisplayRole); - scaleBItem->setData(image_scale_b, ScaleSortRole); - } else { - scaleBItem->setData("N/A", Qt::DisplayRole); - scaleBItem->setData(ScaleNotAvailableSortValue, ScaleSortRole); - } - rowItems.append(scaleBItem); - tableModel->appendRow(rowItems); } @@ -184,10 +172,6 @@ void JFJochViewerImageListWindow::datasetLoaded(std::shared_ptrimage_scale_cc.size() > i) image_scale_cc_percent = dataset->image_scale_cc[i] * 100.0; - double image_scale_b = NAN; - if (dataset->image_scale_b.size() > i) - image_scale_b = dataset->image_scale_b[i]; - addDataRow(i + 1, bkg_estimate, indexing_result, @@ -195,8 +179,7 @@ void JFJochViewerImageListWindow::datasetLoaded(std::shared_ptrUnits("deg"); - data_file.SaveVector("/entry/MX/imageScaleBFactor", image_scale_b_factor.vec())->Units("Angstrom^2"); } } diff --git a/writer/HDF5DataFilePluginMX.h b/writer/HDF5DataFilePluginMX.h index 745c7580..3f0265b1 100644 --- a/writer/HDF5DataFilePluginMX.h +++ b/writer/HDF5DataFilePluginMX.h @@ -61,7 +61,6 @@ class HDF5DataFilePluginMX : public HDF5DataFilePlugin { AutoIncrVector image_scale_factor{NAN}; AutoIncrVector image_scale_cc{NAN}; AutoIncrVector image_scale_mosaicity{NAN}; - AutoIncrVector image_scale_b_factor{NAN}; public: explicit HDF5DataFilePluginMX(const StartMessage& msg); void OpenFile(HDF5File &data_file, const DataMessage& msg, size_t images_per_file) override; diff --git a/writer/HDF5NXmx.cpp b/writer/HDF5NXmx.cpp index c631f167..d0ff03e5 100644 --- a/writer/HDF5NXmx.cpp +++ b/writer/HDF5NXmx.cpp @@ -1037,7 +1037,6 @@ void NXmx::EndResultVectors(const EndMessage &end) { SaveVectorIfMissing(*hdf5_file, "/entry/MX/imageScaleFactor", end.image_scale_factor); SaveVectorIfMissing(*hdf5_file, "/entry/MX/imageScaleCC", end.image_scale_cc); SaveVectorIfMissing(*hdf5_file, "/entry/MX/imageScaleMosaicity", end.image_scale_mosaicity, "deg"); - SaveVectorIfMissing(*hdf5_file, "/entry/MX/imageScaleBFactor", end.image_scale_b_factor, "Angstrom^2"); if (!end.niggli_class.empty()) SaveVectorIfMissing(*hdf5_file, "/entry/MX/niggliClass", end.niggli_class); }