diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 9d1db4f6..f261d670 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -6,6 +6,8 @@ This is an UNSTABLE release. It includes many experimental features, as well as * Space-group search: a screw axis is now claimed from how decisively its predicted-absent reflections are weaker than the rest of their own axial row, instead of from a minimum count of them, so a screw survives a sweep that recorded few axial reflections and is refused on a row too weak to decide either way. * Bragg integration: the profile fit's `background_variance` now takes the fitted intensity itself out of the fit variance instead of `max(0, I)`, so a reflection that fluctuated below zero no longer reports a background variance two to three times too small and is no longer weighted up for it. * Scaling: the rotation merge weights each combined full by its variance rebuilt at the reflection's mean intensity rather than by the full's own sigma, as the stills merge already did. +* rugnux: Rotation data are integrated on **every frame whose spots the sweep's lattice explains**, instead of only on frames that would also index on their own; the reported indexing rate still counts the latter. +* Scaling: a rotation frame too sparse to fit a rocking width of its own now takes the run's median instead of a fixed default. * rugnux: The error-model **a** and **b** are reported in XDS's convention, and `_reflns.jfjoch_diffrn_ISa` now carries the whole-range `1/sqrt(a*b)` that XDS's ISa denotes; the strong-reflection asymptote moves to `_reflns.jfjoch_diffrn_ISa_asymptotic`. **A file written by an earlier version carries the asymptote under the old name.** * Bragg integration: the background ring's outer radius default changes from 10 px to **13 px**, which roughly doubles the pixels behind each background estimate; the signal disk is unchanged. * Scaling: when too few reflections are strong enough to constrain the error model's systematic term **b**, it is now held at zero and **ISa is reported as unmeasured** rather than as the spurious value the fit would otherwise extrapolate. diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 502a8b75..1e6f3b97 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -452,7 +452,11 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg, ewald_dist_cutoff = experiment.GetBraggIntegrationSettings().GetFixedProfileRadius_recipA().value() * 3.0f; float wedge_deg = 0.0f; - float mos_deg = 0.1f; + // The rocking width, when this frame has one. A frame too sparse to fit its own (and, on rotation, + // one integrated purely from the sweep's lattice) predicts with the default below, but must NOT + // report that default onward: RotationScaleMerge averages the reported values in frame order to + // recompute every partiality, and a placeholder entered there is read as a measurement. + std::optional mos_measured; if (experiment.GetGoniometer().has_value()) { // Full oscillation wedge of one frame; BraggPredictionRot halves it to the +/- half-wedge of the @@ -461,8 +465,8 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg, wedge_deg = experiment.GetGoniometer()->GetWedge_deg(); if (msg.mosaicity_deg) { - mos_deg = msg.mosaicity_deg.value(); - mosaicity[msg.number] = mos_deg; + mos_measured = msg.mosaicity_deg.value(); + mosaicity[msg.number] = *mos_measured; } // Second pass of the rotation two-pass: widen the prediction to the frame-order-smoothed mosaicity // that RotationScaleMerge fitted in the first pass. Take the MAX with this frame's own estimate so @@ -471,15 +475,16 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg, if (msg.number >= 0 && msg.number < static_cast(prediction_mosaicity_override_.size()) && std::isfinite(prediction_mosaicity_override_[msg.number]) && prediction_mosaicity_override_[msg.number] > 0.0f) { - mos_deg = std::max(mos_deg, prediction_mosaicity_override_[msg.number]); - mosaicity[msg.number] = mos_deg; + mos_measured = std::max(mos_measured.value_or(0.0f), prediction_mosaicity_override_[msg.number]); + mosaicity[msg.number] = *mos_measured; } } + const float mos_deg = mos_measured.value_or(0.1f); IntegrationOutcome i_outcome{ .geom = outcome.experiment.GetDiffractionGeometry(), .latt = latt, - .mosaicity_deg = mos_deg, + .mosaicity_deg = mos_measured, .image_scale_cc = msg.image_scale_cc, }; @@ -556,6 +561,8 @@ IndexAndRefine::DetermineRefineAnalyze(DataMessage &msg, const SpotFindingSettin if (!outcome.lattice_candidate.has_value()) return std::nullopt; + // AnalyzeIndexing answers "is this frame worth integrating"; msg.indexing_result carries the + // stricter "does this frame index on its own", which on rotation is not the same question. if (!AnalyzeIndexing(msg, outcome.experiment, *outcome.lattice_candidate, outcome.extra_lattice_candidates)) return std::nullopt; @@ -577,7 +584,11 @@ void IndexAndRefine::ProcessImage(DataMessage &msg, } bool IndexAndRefine::IndexFrameOnly(DataMessage &msg, const SpotFindingSettings &spot_finding_settings) { - return DetermineRefineAnalyze(msg, spot_finding_settings).has_value(); + // The per-frame verdict, not "was there anything to integrate": this is what the rotation first pass + // scores candidate lattices on, and on rotation a frame can be integrable without indexing on its + // own (AnalyzeIndexing), which would score a sparse frame for every candidate alike. + return DetermineRefineAnalyze(msg, spot_finding_settings).has_value() + && msg.indexing_result.value_or(false); } std::optional IndexAndRefine::FinalizeRotationIndexing() { diff --git a/image_analysis/indexing/AnalyzeIndexing.cpp b/image_analysis/indexing/AnalyzeIndexing.cpp index 22b10249..056b8fe9 100644 --- a/image_analysis/indexing/AnalyzeIndexing.cpp +++ b/image_analysis/indexing/AnalyzeIndexing.cpp @@ -427,10 +427,22 @@ bool AnalyzeIndexing(DataMessage &message, // once over the sweep and then integrate every frame from that lattice - none of them re-decides // per frame whether a frame may be integrated. constexpr float min_frac = 0.20f; - if (nspots_indexed >= viable_cell_min_spots && nspots_indexed >= std::lround(min_frac * nspots_ref)) { + const bool lattice_fits = nspots_indexed >= std::lround(min_frac * nspots_ref); + // Two different questions. "Does this frame index?" - reported as the indexing rate, and what the + // rotation first pass scores candidate lattices on - needs the absolute floor too, because a + // handful of spots sit on almost any lattice by chance. "Is this frame worth integrating?" needs + // only the consistency part, and only where the lattice does not come from this frame: on rotation + // it is the whole sweep's and the orientation comes from the goniometer, so a sparse frame whose + // few spots all lie on it is a frame of the same crystal, not an unindexed one. Refusing it throws + // away every reflection it records - on a weakly diffracting crystal, most of the dataset, which is + // why XDS's INTEGRATE and dials.integrate predict on every image of the sweep. A frame whose spots + // largely MISS the lattice is a different matter and is still refused. + const bool frame_indexes = lattice_fits && nspots_indexed >= viable_cell_min_spots; + const bool integrate_frame = experiment.IsRotationIndexing() ? lattice_fits : frame_indexes; + if (integrate_frame) { auto uc = latt.GetUnitCell(); if (ok(uc.a) && ok(uc.b) && ok(uc.c) && ok(uc.alpha) && ok(uc.beta) && ok(uc.gamma)) { - message.indexing_result = true; + message.indexing_result = frame_indexes; indexing_lattice_count++; assert(indexed_spots.size() == message.spots.size()); @@ -507,6 +519,6 @@ bool AnalyzeIndexing(DataMessage &message, auto end_time = std::chrono::steady_clock::now(); message.index_analysis_time_s = std::chrono::duration(end_time - start_time).count(); message.indexing_lattice_count = indexing_lattice_count; - message.indexing_result = outcome; + message.indexing_result = outcome && frame_indexes; return outcome; } diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 94439bcd..71276f1d 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -649,6 +649,20 @@ void RotationScaleMerge::SmoothMosaicityAndPartiality() { for (int o = 0; o < n_frames; ++o) mos_smooth[o] = static_cast(mos_raw[o]); } + // A frame too sparse to fit its own rocking width still records reflections, and on a weakly + // diffracting crystal most frames are that sparse - so the smoothing window can be empty over long + // stretches. Those frames used to fall back on the fixed default, a width that has nothing to do + // with this crystal; the run's own median is the same measurement its neighbours made. + std::vector measured; + for (int o = 0; o < n_frames; ++o) + if (std::isfinite(mos_smooth[o])) measured.push_back(mos_smooth[o]); + if (!measured.empty()) { + std::nth_element(measured.begin(), measured.begin() + measured.size() / 2, measured.end()); + const auto median = static_cast(measured[measured.size() / 2]); + for (int o = 0; o < n_frames; ++o) + if (!std::isfinite(mos_smooth[o])) mos_smooth[o] = median; + } + // Recompute each partial's partiality from the smoothed mosaicity (same wedge the predictor used). // Frames without a mosaicity keep the stored partiality. const double wedge = gon ? std::fabs(gon->GetWedge_deg()) : 0.0;