Files
Jungfraujoch/common/ScalingSettings.cpp
leonarski_fandClaude Opus 5 a6be35ccdb Azimuthal integration: optional sigma clipping of the reported profile
The profile is the MEAN of each bin, so a few strong reflections landing in a bin
lift it exactly as a smooth powder ring does. That is the wrong quantity whenever
the profile is wanted as a background rather than as a measurement of what is in
the bin - the ice score being the case in point, where reading a plain profile
INVERTED the metric: over 37 rotation crystals the two highest-scoring crystals
had no ice at all.

The adaptive spot finder already computes the right thing, a sigma-clipped
per-resolution-ring background, as a byproduct of its own threshold. Where it
runs, the ice score uses that. Where it does not - --no-adaptive-spots,
--azint-only, and anything reading the profile the broker wrote - there was no way
to get it. This adds one: azim_int_settings.sigma_clip (rugnux --azim-sigma-clip),
0 = off, minimum 2 because a tighter clip rejects a large part of a clean Gaussian
bin and biases the estimate low rather than removing outliers.

Two clip passes follow the plain one, matching the finder's recipe - the first
pass's standard deviation is itself inflated by the peaks being removed, so one
pass leaves a threshold that is still too generous. A bin with fewer than eight
pixels is left alone: at the detector edge and behind the beam stop there is no
spread to clip on.

Both engines do it. On the GPU the accept range is computed by a small kernel and
stays resident, so a clip pass is one more read of the same pixels and no round
trip; the two accumulation kernels take the range as a pointer that is null on the
plain pass. Measured on a JUNGFRAU rotation dataset, non-adaptive path: azimuthal
integration 0.02 -> 0.06 ms per image, exactly the 3x the extra passes predict,
against a 0.34 ms per-image total.

Note what the result IS: the smooth background under the peaks, not the bin mean.
It should not be switched on where a ring's integrated intensity is wanted - the
powder-ring geometry fit reads ring peaks, and those are what a clip is designed
to remove. Off by default, so nothing changes unless it is asked for.

Not exposed over the REST API - that needs the generated model regenerated, which
is a separate step.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
2026-08-07 00:18:52 +02:00

295 lines
8.7 KiB
C++

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "ScalingSettings.h"
ScalingSettings& ScalingSettings::MergeFriedel(bool input) {
merge_friedel = input;
return *this;
}
ScalingSettings& ScalingSettings::HighResolutionLimit_A(double limit) {
if (limit <= 0.0)
throw JFJochException(JFJochExceptionCategory::InputParameterBelowMin, "High resolution limit must be positive");
high_resolution_limit_A = limit;
return *this;
}
ScalingSettings& ScalingSettings::HighResolutionLimit_A(std::optional<double> limit) {
if (limit.has_value() && limit.value() <= 0.0)
throw JFJochException(JFJochExceptionCategory::InputParameterBelowMin, "High resolution limit must be positive");
high_resolution_limit_A = limit;
return *this;
}
bool ScalingSettings::GetMergeFriedel() const {
return merge_friedel;
}
ScalingSettings &ScalingSettings::RefineRotationWedge(bool input) {
refine_wedge = input;
return *this;
}
bool ScalingSettings::GetRefineWedge() const {
return refine_wedge;
}
std::optional<double> ScalingSettings::GetHighResolutionLimit_A() const {
return high_resolution_limit_A;
}
double ScalingSettings::GetMinMosaicity() const {
return 0.001;
}
double ScalingSettings::GetMaxMosaicity() const {
return 1.0;
}
double ScalingSettings::GetMinWedge() const {
return 0.001;
}
double ScalingSettings::GetMaxWedge() const {
return 10.0;
}
double ScalingSettings::GetDefaultMosaicity() const {
return 0.1;
}
ScalingSettings &ScalingSettings::RotationWedgeForScaling(std::optional<double> input) {
if (input) {
// TODO: Use fmt
if (input.value() < GetMinWedge() || input.value() > GetMaxWedge())
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Wedge for scaling must be between " + std::to_string(GetMinWedge()) +
" and " + std::to_string(GetMaxWedge()));
}
wedge_for_scaling = input;
return *this;
}
std::optional<double> ScalingSettings::GetRotationWedgeForScaling() const {
return wedge_for_scaling;
}
ScalingSettings &ScalingSettings::MinPartiality(double input) {
if (min_partiality < 0.0 || min_partiality > 1.0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Min partiality must be between 0 and 1");
min_partiality = input;
return *this;
}
double ScalingSettings::GetMinCCForImage() const {
return min_cc_for_image;
}
double ScalingSettings::GetSearchMinZeta() const {
return search_min_zeta;
}
ScalingSettings &ScalingSettings::SearchMinZeta(double input) {
if (input < 0.0 || input >= 1.0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Search zeta limit must be in [0,1)");
search_min_zeta = input;
return *this;
}
ScalingSettings &ScalingSettings::MinCCForImage(double input) {
if (input < 0.0 || input > 1.0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Min CC for image must be between 0 and 1");
min_cc_for_image = input;
return *this;
}
double ScalingSettings::GetOutlierRejectNsigma() const {
return outlier_reject_nsigma;
}
ScalingSettings &ScalingSettings::OutlierRejectNsigma(double input) {
outlier_reject_nsigma = input; // <= 0 disables; no upper bound (large = effectively off)
return *this;
}
ScalingSettings &ScalingSettings::ScaleFulls(bool input) {
scale_fulls = input;
return *this;
}
bool ScalingSettings::GetScaleFulls() const {
return scale_fulls;
}
ScalingSettings &ScalingSettings::AbsorptionIter(int input) {
if (input < 0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Absorption iterations must be non-negative");
absorption_iter = input;
return *this;
}
int ScalingSettings::GetAbsorptionIter() const {
return absorption_iter;
}
ScalingSettings &ScalingSettings::CorrectionSurfaces(bool input) {
correction_surfaces = input;
return *this;
}
bool ScalingSettings::GetCorrectionSurfaces() const {
return correction_surfaces;
}
ScalingSettings &ScalingSettings::StillsPartialityRefine(bool input) {
stills_partiality_refine = input;
return *this;
}
bool ScalingSettings::GetStillsPartialityRefine() const {
return stills_partiality_refine;
}
ScalingSettings &ScalingSettings::ExpectedVarianceMerge(bool input) {
expected_variance_merge = input;
return *this;
}
bool ScalingSettings::GetExpectedVarianceMerge() const {
return expected_variance_merge;
}
float ScalingSettings::GetIceMinScore() const {
return ice_min_score;
}
float ScalingSettings::GetIceMinSpotRatio() const {
return ice_min_spot_ratio;
}
ScalingSettings &ScalingSettings::IceMinSpotRatio(float input) {
if (input < 0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Ice spot-ratio gate must be non-negative");
ice_min_spot_ratio = input;
return *this;
}
ScalingSettings &ScalingSettings::IceMinScore(float input) {
if (input < 0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Ice score gate must be non-negative");
ice_min_score = input;
return *this;
}
ScalingSettings &ScalingSettings::SmoothGDegrees(double input) {
if (input < 0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Smooth-G range must be non-negative");
smooth_g_deg = input;
return *this;
}
double ScalingSettings::GetSmoothGDegrees() const {
return smooth_g_deg;
}
ScalingSettings &ScalingSettings::RelativeBDegrees(double input) {
if (input < 0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Relative-B batch width must be non-negative");
relative_b_deg = input;
return *this;
}
double ScalingSettings::GetRelativeBDegrees() const {
return relative_b_deg;
}
double ScalingSettings::GetMinPartiality() const {
return min_partiality;
}
ScalingSettings &ScalingSettings::ForcedMosaicity(std::optional<double> input) {
if (input.has_value() && (input.value() < GetMinMosaicity() || input.value() > GetMaxMosaicity()))
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Forced mosaicity must be between " + std::to_string(GetMinMosaicity()) +
" and " + std::to_string(GetMaxMosaicity()));
forced_mosaicity = input;
return *this;
}
std::optional<double> ScalingSettings::GetForcedMosaicity() const {
return forced_mosaicity;
}
ScalingSettings &ScalingSettings::CaptureUncertaintyCoeff(double input) {
if (input < 0.0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Capture uncertainty coefficient must be non-negative");
capture_uncertainty_coeff = input;
return *this;
}
double ScalingSettings::GetCaptureUncertaintyCoeff() const {
return capture_uncertainty_coeff;
}
ScalingSettings &ScalingSettings::MinCapturedFraction(double input) {
if (input < 0.0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Minimum captured fraction must be non-negative");
min_captured_fraction = input;
return *this;
}
double ScalingSettings::GetMinCapturedFraction() const {
return min_captured_fraction;
}
ScalingSettings &ScalingSettings::RfreeFraction(double input) {
if (input < 0.0 || input > 1.0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "R-free fraction must be between 0 and 1");
rfree_fraction = input;
return *this;
}
double ScalingSettings::GetRfreeFraction() const {
return rfree_fraction;
}
ScalingSettings &ScalingSettings::ResolutionCutoff(ResolutionCutoffMethod input) {
resolution_cutoff = input;
return *this;
}
ResolutionCutoffMethod ScalingSettings::GetResolutionCutoff() const {
return resolution_cutoff;
}
ScalingSettings &ScalingSettings::ResolutionCCTarget(double input) {
if (input <= 0.0 || input >= 1.0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Resolution CC target must be between 0 and 1");
resolution_cc_target = input;
return *this;
}
double ScalingSettings::GetResolutionCCTarget() const {
return resolution_cc_target;
}
ScalingSettings &ScalingSettings::ReportShellCount(int input) {
if (input < 1)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Number of report shells must be at least 1");
report_shell_count = input;
return *this;
}
int ScalingSettings::GetReportShellCount() const {
return report_shell_count;
}