Files
Jungfraujoch/common/IndexingSettings.cpp
T
leonarski_fandClaude Opus 4.8 9415c0483b Add opt-in local-SNR spot gate and acceptance-fraction knob (serial stills)
Two opt-in tools for weak serial-stills tuning; both default-off, so the
default pipeline is bit-identical (verified: NmHR 7.85% identical to HEAD).

--local-snr <sigma> (AdaptiveSpotFinderCPU::FilterByLocalSNR): after the loose
per-ring adaptive threshold builds connected-component spots, drop any spot that
does not stand this many sigmas above its OWN LOCAL background (robust median/MAD
of a square annulus), not just the azimuthal ring mean. On structured-background
(XFEL) frames the ring mean underestimates the local diffuse level in some
sectors, so the ring threshold floods; a real Bragg peak still stands many local
sigmas proud. Validated on XFEL stills to separate real peaks from flood at the
pixel level (real median local-SNR ~70 vs flood ~2.6; SNR>=5 keeps ~99.8% of
real peaks, ~14% of flood). GPU-portable (a per-spot local reduction). NOTE: on
the current battery (NmHR/OCP) it is index-rate/CC1/2 neutral -- the flood that
survives as CC clusters overlaps weak-real spots, and only lattice-fit separates
those -- but it is the correct tool for genuinely floody data (ice/jet/loosened
detector) and the right substrate for the online FPGA path.

--min-indexed-fraction <f>: exposes the previously hardcoded 0.20 minimum
indexed-spot fraction (AnalyzeIndexing) as a per-run setting. Lowering it admits
weaker/sparser crystals; on flooded XFEL data the extra lattices are spurious
(pair with --min-image-cc to gate them), on clean synchrotron data there are no
marginal frames so it is a no-op -- useful as a gating-experiment primitive.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
2026-07-24 11:12:22 +02:00

236 lines
7.2 KiB
C++

// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "IndexingSettings.h"
#include "JFJochException.h"
#include "CUDAWrapper.h"
#include <cmath>
#define check_max(param, val, max) if ((val) > (max)) throw JFJochException(JFJochExceptionCategory::InputParameterAboveMax, param)
#define check_min(param, val, min) if ((val) < (min)) throw JFJochException(JFJochExceptionCategory::InputParameterBelowMin, param)
#define check_finite(param, val) if (!std::isfinite(val)) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, param)
IndexingSettings::IndexingSettings() {
if (get_gpu_count() > 0)
algorithm = IndexingAlgorithmEnum::FFBIDX;
else
algorithm = IndexingAlgorithmEnum::None;
}
IndexingSettings &IndexingSettings::ViableCellMinSpots(int64_t input) {
check_min("ViableCellMinSpots", input, 6);
viable_cell_min_spots = input;
return *this;
}
int64_t IndexingSettings::GetViableCellMinSpots() const {
return viable_cell_min_spots;
}
IndexingSettings &IndexingSettings::MinIndexedSpotFraction(float input) {
min_indexed_spot_fraction = input;
return *this;
}
float IndexingSettings::GetMinIndexedSpotFraction() const {
return min_indexed_spot_fraction;
}
IndexingSettings &IndexingSettings::Algorithm(IndexingAlgorithmEnum input) {
switch (input) {
case IndexingAlgorithmEnum::Auto:
case IndexingAlgorithmEnum::FFBIDX:
case IndexingAlgorithmEnum::FFT:
case IndexingAlgorithmEnum::FFTW:
case IndexingAlgorithmEnum::None:
algorithm = input;
break;
default:
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Invalid value for indexing algorithm enum parameter");
}
return *this;
}
IndexingSettings &IndexingSettings::FFT_MaxUnitCell_A(float input) {
check_finite("FFT indexing max unit cell (A)", input);
check_min("FFT indexing max unit cell (A)", input, 50);
check_max("FFT indexing max unit cell (A)", input, 500);
fft_max_unit_cell_A = input;
return *this;
}
IndexingSettings &IndexingSettings::FFT_MinUnitCell_A(float input) {
check_finite("FFT indexing min unit cell (A)", input);
check_min("FFT indexing min unit cell (A)", input, 5);
check_max("FFT indexing min unit cell (A)", input, 40);
fft_min_unit_cell_A = input;
return *this;
}
IndexingSettings & IndexingSettings::FFT_MaxAngle_deg(float input) {
check_finite("FFT indexing max angle (deg)", input);
check_min("FFT indexing max angle (deg)", input, 0);
check_max("FFT indexing max angle (deg)", input, 180);
fft_max_angle_deg = input;
return *this;
}
IndexingSettings & IndexingSettings::FFT_MinAngle_deg(float input) {
check_finite("FFT indexing min angle (deg)", input);
check_min("FFT indexing min angle (deg)", input, 0);
check_max("FFT indexing min angle (deg)", input, 180);
fft_min_angle_deg = input;
return *this;
}
float IndexingSettings::GetFFT_MinAngle_deg() const {
return fft_min_angle_deg;
}
float IndexingSettings::GetFFT_MaxAngle_deg() const {
return fft_max_angle_deg;
}
IndexingSettings &IndexingSettings::FFT_NumVectors(int64_t input) {
check_min("FFT indexing number of search vectors", input, 128);
fft_num_vectors = input;
return *this;
}
IndexingSettings &IndexingSettings::FFT_HighResolution_A(float input) {
check_finite("FFT indexing high resolution (A)", input);
check_min("FFT indexing high resolution (A)", input, 0.5);
check_max("FFT indexing high resolution (A)", input, 6.0);
fft_high_resolution_A = input;
return *this;
}
IndexingAlgorithmEnum IndexingSettings::GetAlgorithm() const {
return algorithm;
}
float IndexingSettings::GetFFT_MaxUnitCell_A() const {
return fft_max_unit_cell_A;
}
float IndexingSettings::GetFFT_MinUnitCell_A() const {
return fft_min_unit_cell_A;
}
int64_t IndexingSettings::GetFFT_NumVectors() const {
return fft_num_vectors;
}
float IndexingSettings::GetFFT_HighResolution_A() const {
return fft_high_resolution_A;
}
IndexingSettings &IndexingSettings::Tolerance(float input) {
check_min("Indexing tolerance", input, 0.0);
check_max("Indexing tolerance", input, 0.5);
indexing_tolerance = input;
return *this;
}
float IndexingSettings::GetTolerance() const {
return indexing_tolerance;
}
int64_t IndexingSettings::GetIndexingThreads() const {
return indexing_threads;
}
IndexingSettings &IndexingSettings::IndexingThreads(int64_t input) {
check_min("Indexing thread count", input, 1);
check_max("Indexing thread count", input, 64);
indexing_threads = input;
return *this;
}
IndexingSettings &IndexingSettings::UnitCellDistTolerance(float input) {
check_min("Relative unit cell distance tolerance vs. reference", input, 0.0001);
check_max("Relative unit cell distance tolerance vs. reference", input, 0.2001);
unit_cell_dist_tolerance_vs_reference = input;
return *this;
}
float IndexingSettings::GetUnitCellDistTolerance() const {
return unit_cell_dist_tolerance_vs_reference;
}
float IndexingSettings::GetUnitCellAngleTolerance_deg() const {
return unit_cell_angle_tolerance_deg;
}
GeomRefinementAlgorithmEnum IndexingSettings::GetGeomRefinementAlgorithm() const {
return refinement;
}
IndexingSettings &IndexingSettings::GeomRefinementAlgorithm(GeomRefinementAlgorithmEnum input) {
refinement = input;
return *this;
}
IndexingSettings & IndexingSettings::IndexIceRings(bool input) {
index_ice_rings = input;
return *this;
}
IndexingSettings & IndexingSettings::RotationIndexing(bool input) {
enable_rotation_indexing = input;
return*this;
}
IndexingSettings & IndexingSettings::RotationIndexingMinAngularRange_deg(float input) {
check_finite("Rotation indexing minimum angular range (deg.)", input);
check_min("Rotation indexing minimum angular range (deg.)", input, 1.0);
rotation_indexing_min_angular_range_deg = input;
return *this;
}
IndexingSettings & IndexingSettings::RotationIndexingAngularStride_deg(float input) {
check_finite("Rotation indexing angular stride (deg.)", input);
check_min("Rotation indexing angular stride (deg.)", input, 0.0);
rotation_indexing_angular_stride_deg = input;
return *this;
}
bool IndexingSettings::GetRotationIndexing() const {
return enable_rotation_indexing;
}
float IndexingSettings::GetRotationIndexingMinAngularRange_deg() const {
return rotation_indexing_min_angular_range_deg;
}
float IndexingSettings::GetRotationIndexingAngularStride_deg() const {
return rotation_indexing_angular_stride_deg;
}
bool IndexingSettings::GetIndexIceRings() const {
return index_ice_rings;
}
bool IndexingSettings::GetBlockingBehavior() const {
return blocking_behavior;
}
IndexingSettings &IndexingSettings::BlockingBehavior(bool input) {
blocking_behavior = input;
return *this;
}
int64_t IndexingSettings::GetMaxExtraLattices() const {
return max_extra_lattices;
}
IndexingSettings &IndexingSettings::MaxExtraLattices(int64_t input) {
check_min("Max extra lattices", input, 0);
check_max("Max extra lattices", input, 10);
max_extra_lattices = input;
return *this;
}