From 4895dc10187ffe6a08748dc0fdbbfc1b1cfd1a8c Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 28 Jul 2026 11:45:49 +0200 Subject: [PATCH] Rotation: let --min-image-cc drop frames that disagree with the merged reference The flag was accepted on rotation data and did nothing - it is read only by the stills merge (Merge.cpp), and the CLI warned about that rather than fixing it. Meanwhile RotationScaleMerge already COMPUTES a per-frame correlation against the merged reference and writes it to the per-image table; nothing acted on it. Wire the two together. A rejected frame has its partials' corr set to 0, which is how a frame already leaves the pipeline - every consumer requires corr > 0, so the combine, the merge and the error model all drop it together. The GPU path reuses the SmoothCorr kernel with a ratio of 0, so one implementation covers both. Off by default (0), and verified bit-identical to the previous binary when off. What it catches, on the two rotation datasets that have a population to catch: a two-lattice crystal - two lattices in two physical AREAS of the sample, so the sweep passes from one to the other and whole blocks of frames measure a different crystal from the one being merged (frames 500-700 index perfectly well at a per-frame CC of 0.22 against 0.47-0.56 either side, in 11 contiguous runs). R_meas 28.6 -> 24.6%, CC1/2 93.6 -> 95.1, high-shell CC 23.4 -> 38.3. a second dataset with 9.5% of frames below CC 0.30: R_meas 24.3 -> 23.4%, CC1/2 92.6 -> 93.4. The criterion is "this frame disagrees with the merged reference", NOT "this frame is off-crystal". It happens to catch both, because a frame that measures nothing and a frame that measures a DIFFERENT crystal fail the same test, and it does not need to know which. For the two-area case that is a workaround, not a treatment: it recovers one crystal by discarding the other, where processing the two as separate sweeps would keep both. The frame-block structure is clean enough that such a split could be detected automatically. WHY THERE IS NO DEFAULT. The per-frame CC is not comparable between datasets - it is as much a measure of data quality as of frame validity. Measured medians across the battery run from 0.30 to 0.81, so one absolute bound removes 13 frames from one dataset and 584 of 1800 from another: battery at --min-image-cc 30, 33 crystals: no point group changed (30/33), four crystals clearly better (one +5.4 CC1/2 points, the two-lattice case above, and ISa gains of 1.3-4.6 on three others) - and one healthy crystal lost a third of its frames and with them its high-resolution shell (CC1/2_hi 26.2 -> 2.0). This is the same trap as an absolute bound on any per-operator or per-frame agreement statistic, and the same one the per-frame scale guard avoids by measuring against the run's own median. A principled version would cut on the SHAPE of the per-frame CC distribution - a dataset with a bad subpopulation is bimodal, a uniformly weak one is not - rather than on an absolute value. Until that exists this stays opt-in, and the per-image CC it keys on is already in the _image.dat table for anyone choosing a value. Co-Authored-By: Claude Opus 5 (1M context) --- .../scale_merge/RotationScaleMerge.cpp | 33 ++++++++++++++++++- .../scale_merge/RotationScaleMerge.h | 6 ++++ rugnux/rugnux_cli.cpp | 4 --- 3 files changed, 38 insertions(+), 5 deletions(-) diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index b37f3237..c8e3f4ae 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -29,7 +29,8 @@ namespace { // These mirror the per-image ScaleOnTheFly / Merge rocking-curve physics verbatim so this flat // implementation is numerically identical - see the comments there for the details. - constexpr size_t MIN_REFLECTIONS = 20; // per-frame scale needs at least this many + constexpr size_t MIN_REFLECTIONS = 20; + constexpr int64_t MIN_REFLECTIONS_FOR_IMAGE_CC = 20; // below this a frame's CC means nothing // per-frame scale needs at least this many constexpr double SCALE_ROBUST_K = 3.0; // Cauchy loss scale (sigma units) for the per-frame G fit // A fitted per-frame scale below this fraction of the run's median is not a measurement of anything: @@ -175,6 +176,7 @@ RotationScaleMerge::RotationScaleMerge(const DiffractionExperiment &experiment, merge_friedel = s.GetMergeFriedel(); capture_uncertainty_coeff = s.GetCaptureUncertaintyCoeff(); min_captured_fraction = s.GetMinCapturedFraction(); + min_cc_for_image = s.GetMinCCForImage(); reject_nsigma = s.GetOutlierRejectNsigma(); reject_outliers = reject_nsigma > 0.0; rfree_fraction = s.GetRfreeFraction(); @@ -1993,6 +1995,35 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, } FinalizePerFrameScale(cc, cc_n, partial_scaled); + // --- 2b. Drop frames that do not agree with the merged reference (--min-image-cc). --- + if (min_cc_for_image > 0.0) { + std::vector reject(n_frames, 0); + int n_rejected = 0; + for (int f = 0; f < n_frames; ++f) + if (std::isfinite(cc[f]) && cc_n[f] >= MIN_REFLECTIONS_FOR_IMAGE_CC && cc[f] < min_cc_for_image) { + reject[f] = 1; + ++n_rejected; + } + if (n_rejected > 0) { + // Zeroing corr is how a frame leaves the pipeline: every consumer requires corr > 0, so the + // frame's partials stop being usable for the combine, the merge and the error model alike. + bool rejected_on_gpu = false; +#ifdef JFJOCH_USE_CUDA + if (gpu_active_) { + std::vector ratio(n_frames, 1.0); + for (int f = 0; f < n_frames; ++f) if (reject[f]) ratio[f] = 0.0; + gpu_->SmoothCorr(reject.data(), ratio.data()); + rejected_on_gpu = true; + } +#endif + if (!rejected_on_gpu) + for (auto &o : partials) + if (reject[o.frame]) o.corr = 0.0f; + logger.Info("Rejected {} of {} frames correlating below {:.2f} with the merged reference", + n_rejected, n_frames, min_cc_for_image); + } + } + // --- 3. 3D combine of per-frame partials into fulls (fulls inherit their ASU group here). --- bool combined_on_gpu = false; bool scaled_fulls_on_gpu = false; diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 2dd3be38..adb37fa9 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -100,6 +100,12 @@ private: bool merge_friedel = true; double capture_uncertainty_coeff = 0.0; double min_captured_fraction = 0.0; + + // Drop a frame's observations entirely when the frame disagrees with the merged reference below this + // correlation (--min-image-cc). A mis-centred or off-crystal frame still produces spots, still + // indexes and still integrates - it just measures something that is not the crystal's diffraction, + // and nothing downstream removes it. 0 = off. + double min_cc_for_image = 0.0; double reject_nsigma = 0.0; bool reject_outliers = false; double rfree_fraction = 0.0; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 383625d9..d1e1e66f 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -1425,10 +1425,6 @@ int main(int argc, char **argv) { scaling_settings.CaptureUncertaintyCoeff(capture_uncertainty_arg.value_or(rotation_indexing ? 1.0 : 0.0)); scaling_settings.ForcedMosaicity(forced_mosaicity_arg); scaling_settings.MinCCForImage(min_image_cc / 100.0); // --min-image-cc is in percent; the setting is a fraction - // The per-image CC filter is consumed by the stills merge (MergeOnTheFly) only - RotationScaleMerge - // never reads it - so on rotation data it would otherwise be accepted and silently do nothing. - if (min_image_cc > 0.0 && rotation_indexing) - logger.Warning("--min-image-cc is ignored for rotation data: it filters the stills merge only"); scaling_settings.OutlierRejectNsigma( outlier_reject_nsigma.value_or(rotation_indexing ? REJECT_OUTLIERS_DEFAULT_NSIGMA : 0.0));