diff --git a/rugnux/RugnuxFirstPass.cpp b/rugnux/RugnuxFirstPass.cpp index 599a59494..f3423f432 100644 --- a/rugnux/RugnuxFirstPass.cpp +++ b/rugnux/RugnuxFirstPass.cpp @@ -1436,6 +1436,96 @@ bool Rugnux::FirstPassRotationIndexing(PipelineLocals &p) { } } + // Long-axis rescue. When the de-novo cell indexes few validation frames, a long, finely-spaced + // axis was likely lost: the unconstrained FFT either collapsed it to a short sub-multiple or let + // a denser supercell over-fit the accumulated cloud (a small global-orientation error throws the + // many high-order reflections off along the fine axis, so the true cell scores worst on the raw + // cloud). Recover the true metric with a second search whose FFT unit-cell bound is widened past + // the default, then RE-INDEX at full resolution constrained by the cell that comes back, as a + // reference (the -C path): the reference filter drops the collapsed/supercell candidates and + // refines an accurate global lattice. Only runs after a poor standard pass, so well-indexing + // crystals never pay for it and their result is untouched. Returns the re-index where it indexes + // more validation frames than the `standing` pass it was asked about. + auto max_axis = [](const RotationIndexerResult &r) { + const auto uc = r.lattice.GetUnitCell(); + return std::max({uc.a, uc.b, uc.c}); + }; + const auto long_axis_rescue = [&](const FirstPass &standing) -> std::optional { + if (cancelled_ || !standing.result.has_value() + || standing.score >= 0.5 * static_cast(validation.size())) + return std::nullopt; + auto coarse_settings = experiment_.GetIndexingSettings(); + // The resolution setting only sizes the FFT histogram - it does not filter which spots go + // into it - so this is here to pay for the widened bound below, not to make the long axis + // easier to pick out. Measured, coarsening does the opposite: the acceptance angle widens + // (0.159 -> 0.579 deg) but the true axis falls from rank 1996 to 16361 of 16384 directions. + coarse_settings.FFT_HighResolution_A(3.5f); + // The widened bound is what the rescue turns on, and it costs nothing: the histogram length is + // 4*pi*max_length/high_resolution, so coarsening the resolution to 3.5 A has already + // SHRUNK the transform by 1.75x - spending that back on reach leaves the coarse pass at + // 1.14x the standard one. Without it the rescue cannot recover what it exists to recover: + // a cell whose long axis is past the default 500 A bound is unrepresentable, so the + // "recovered" metric would still be a sub-multiple. Only reached after a poor pass, so a + // well-indexing crystal never pays. Stills never reach here. + constexpr float RESCUE_MAX_UNIT_CELL_A = 1000.0f; + if (coarse_settings.GetFFT_MaxUnitCell_A() < RESCUE_MAX_UNIT_CELL_A) + coarse_settings.FFT_MaxUnitCell_A(RESCUE_MAX_UNIT_CELL_A); + coarse_settings.IndexingThreads(2); + coarse_settings.RefineThreads(std::clamp(config_.nthreads / 2, 1, MAX_REFINE_THREADS)); + IndexerThreadPool coarse_pool(coarse_settings, IndexerConstruction::OnFirstUse); + // Coarse first pass: keep the recovered cell with the LONGEST axis directly. Its full- + // resolution per-frame validation would be low (the coarse cell is metrically right but + // imprecise), so do NOT score it here - it is only the reference for the constrained re-index. + std::vector> cris; + cris.reserve(schemes.size()); + for (const auto &[name, ordinals] : schemes) { + if (cancelled_) break; + cris.push_back(feed_scheme(coarse_pool, ordinals)); + } + std::vector> cfut; + cfut.reserve(cris.size()); + for (auto &ri : cris) { + RotationIndexer *rp = ri.get(); + cfut.push_back(std::async(std::launch::async, [rp] { rp->RunIndexing(); })); + } + for (auto &f : cfut) f.get(); + std::optional coarse_ref; + for (auto &ri : cris) { + auto r = ri->GetLattice(); + if (r.has_value() && (!coarse_ref.has_value() || max_axis(*r) > max_axis(*coarse_ref))) + coarse_ref = r; + } + // A recovered axis at least 30% longer than the standard pass found signals a rescued metric. + if (!cancelled_ && coarse_ref.has_value() + && max_axis(*coarse_ref) > 1.3 * max_axis(*standing.result)) { + logger.Info("Long-axis rescue: coarse pass recovered a {:.0f} A axis (was {:.0f} A); " + "re-indexing constrained by that cell", max_axis(*coarse_ref), max_axis(*standing.result)); + experiment_.SetUnitCell(coarse_ref->lattice.GetUnitCell()); + // The re-index must be able to REPRESENT the cell the coarse pass just recovered. The + // standard pool was sized at the default bound, so re-using it here would search for a + // long axis it cannot hold and hand back the sub-multiple the rescue exists to replace. + // count_indexed forces the lattice and never touches the pool, so the existing + // IndexAndRefine is reused; only the FFT needs the wider one. + auto wide_settings = experiment_.GetIndexingSettings(); + wide_settings.IndexingThreads(2); + wide_settings.RefineThreads(std::clamp(config_.nthreads / 2, 1, MAX_REFINE_THREADS)); + RaiseFFTBoundForKnownCell(wide_settings, experiment_); + FirstPass constrained; + if (wide_settings.GetFFT_MaxUnitCell_A() > indexing_bound_A) { + IndexerThreadPool wide_pool(wide_settings, IndexerConstruction::OnFirstUse); + constrained = pick_best(wide_pool, *indexer); + } else { + constrained = pick_best(*indexer_pool, *indexer); + } + experiment_.SetUnitCell(std::nullopt); // leave the space-group / cell determination de-novo + if (constrained.result.has_value() && constrained.score > standing.score) + return constrained; + } + return std::nullopt; + }; + // Set where the beam-centre check below has already asked the rescue at the centre it keeps. + bool long_axis_asked = false; + // The background centre, tried as a SECOND FIRST PASS. // // The estimate itself is free - it is read off the projection --detect-beam-stop has already @@ -1654,20 +1744,37 @@ bool Rugnux::FirstPassRotationIndexing(PipelineLocals &p) { // and then only the pooled excess over chance - each lattice measured against its // own wrong-spindle null, so a larger cell is not credited for fitting more by // chance - says which centre sees it. + // + // Each centre is judged on the best lattice it gives, so the long-axis rescue + // is asked at BOTH before they are compared. The lattices standing here index + // no majority anywhere, which is exactly the pass that rescue exists for, and + // the two counts of spots on two lattices that are not the crystal's are close + // to a coin toss: measured on a 545 A axis, the measured centre (3.6 px out + // across the spindle) won it on 16.7% against 15.5%, the rescue then ran + // there alone and could not see the long axis, and the run shipped a P1 + // sub-cell - while at the file's centre the same rescue indexes 60/60. + if (auto rescued = long_axis_rescue(alt)) + alt = *rescued; PooledEvidence at_measured, at_header; if (alt.result.has_value()) at_measured = pooled_evidence(*indexer, *alt.result); restore_beam_center(header_x, header_y); + if (auto rescued = long_axis_rescue(best)) + best = *rescued; if (best.result.has_value()) at_header = pooled_evidence(*indexer, *best.result); + long_axis_asked = true; if (alt.result.has_value() && beats_chance(at_measured) && ValidationEvidencePrefers(at_header, at_measured)) { try_beam_center(measured_x, measured_y); logger.Warning("Beam centre from the background: neither centre indexes a " - "validation frame on its own, but the pooled spots put " + "majority of the validation frames ({}/{} at the measured, " + "{}/{} at the file's), but the pooled spots put " "{:.1f}% on the lattice at the measured ({:.2f},{:.2f}) against " "{:.1f}% at the file's (chance {:.1f}% and {:.1f}%) - the run " "adopts the measured centre", + alt.score, static_cast(validation.size()), + best.score, static_cast(validation.size()), 100.0 * at_measured.on_lattice / std::max(1, at_measured.spots), measured_x, measured_y, 100.0 * at_header.on_lattice / std::max(1, at_header.spots), @@ -1675,9 +1782,10 @@ bool Rugnux::FirstPassRotationIndexing(PipelineLocals &p) { 100.0 * at_header.by_chance / std::max(1, at_header.spots)); best = alt; } else { - logger.Info("Beam centre check: neither the file's centre ({}/{} frames) nor " - "the measured one {:.2f} px away ({}/{}) indexes this run - the " - "file's centre is kept and this is not a beam-centre failure", + logger.Info("Beam centre check: the file's centre ({}/{} frames) and the " + "measured one {:.2f} px away ({}/{}) - the measured centre puts no " + "more of the pooled spots on its lattice, so the file's centre is " + "kept", best.score, static_cast(validation.size()), moved, alt.score, static_cast(validation.size())); } @@ -1857,89 +1965,10 @@ bool Rugnux::FirstPassRotationIndexing(PipelineLocals &p) { } } - // Long-axis rescue. When the de-novo cell indexes few validation frames, a long, finely-spaced - // axis was likely lost: the unconstrained FFT either collapsed it to a short sub-multiple or let - // a denser supercell over-fit the accumulated cloud (a small global-orientation error throws the - // many high-order reflections off along the fine axis, so the true cell scores worst on the raw - // cloud). Recover the true metric with a second search whose FFT unit-cell bound is widened past - // the default, then RE-INDEX at full resolution constrained by the cell that comes back, as a - // reference (the -C path): the reference filter drops the collapsed/supercell candidates and - // refines an accurate global lattice. Only runs after a poor standard pass, so well-indexing - // crystals never pay for it and their result is untouched. - auto max_axis = [](const RotationIndexerResult &r) { - const auto uc = r.lattice.GetUnitCell(); - return std::max({uc.a, uc.b, uc.c}); - }; - if (!cancelled_ && best.result.has_value() - && best.score < 0.5 * static_cast(validation.size())) { - auto coarse_settings = experiment_.GetIndexingSettings(); - // The resolution setting only sizes the FFT histogram - it does not filter which spots go - // into it - so this is here to pay for the widened bound below, not to make the long axis - // easier to pick out. Measured, coarsening does the opposite: the acceptance angle widens - // (0.159 -> 0.579 deg) but the true axis falls from rank 1996 to 16361 of 16384 directions. - coarse_settings.FFT_HighResolution_A(3.5f); - // The widened bound is what the rescue turns on, and it costs nothing: the histogram length is - // 4*pi*max_length/high_resolution, so coarsening the resolution to 3.5 A has already - // SHRUNK the transform by 1.75x - spending that back on reach leaves the coarse pass at - // 1.14x the standard one. Without it the rescue cannot recover what it exists to recover: - // a cell whose long axis is past the default 500 A bound is unrepresentable, so the - // "recovered" metric would still be a sub-multiple. Only reached after a poor pass, so a - // well-indexing crystal never pays. Stills never reach here. - constexpr float RESCUE_MAX_UNIT_CELL_A = 1000.0f; - if (coarse_settings.GetFFT_MaxUnitCell_A() < RESCUE_MAX_UNIT_CELL_A) - coarse_settings.FFT_MaxUnitCell_A(RESCUE_MAX_UNIT_CELL_A); - coarse_settings.IndexingThreads(2); - coarse_settings.RefineThreads(std::clamp(config_.nthreads / 2, 1, MAX_REFINE_THREADS)); - IndexerThreadPool coarse_pool(coarse_settings, IndexerConstruction::OnFirstUse); - // Coarse first pass: keep the recovered cell with the LONGEST axis directly. Its full- - // resolution per-frame validation would be low (the coarse cell is metrically right but - // imprecise), so do NOT score it here - it is only the reference for the constrained re-index. - std::vector> cris; - cris.reserve(schemes.size()); - for (const auto &[name, ordinals] : schemes) { - if (cancelled_) break; - cris.push_back(feed_scheme(coarse_pool, ordinals)); - } - std::vector> cfut; - cfut.reserve(cris.size()); - for (auto &ri : cris) { - RotationIndexer *rp = ri.get(); - cfut.push_back(std::async(std::launch::async, [rp] { rp->RunIndexing(); })); - } - for (auto &f : cfut) f.get(); - std::optional coarse_ref; - for (auto &ri : cris) { - auto r = ri->GetLattice(); - if (r.has_value() && (!coarse_ref.has_value() || max_axis(*r) > max_axis(*coarse_ref))) - coarse_ref = r; - } - // A recovered axis at least 30% longer than the standard pass found signals a rescued metric. - if (!cancelled_ && coarse_ref.has_value() - && max_axis(*coarse_ref) > 1.3 * max_axis(*best.result)) { - logger.Info("Long-axis rescue: coarse pass recovered a {:.0f} A axis (was {:.0f} A); " - "re-indexing constrained by that cell", max_axis(*coarse_ref), max_axis(*best.result)); - experiment_.SetUnitCell(coarse_ref->lattice.GetUnitCell()); - // The re-index must be able to REPRESENT the cell the coarse pass just recovered. The - // standard pool was sized at the default bound, so re-using it here would search for a - // long axis it cannot hold and hand back the sub-multiple the rescue exists to replace. - // count_indexed forces the lattice and never touches the pool, so the existing - // IndexAndRefine is reused; only the FFT needs the wider one. - auto wide_settings = experiment_.GetIndexingSettings(); - wide_settings.IndexingThreads(2); - wide_settings.RefineThreads(std::clamp(config_.nthreads / 2, 1, MAX_REFINE_THREADS)); - RaiseFFTBoundForKnownCell(wide_settings, experiment_); - FirstPass constrained; - if (wide_settings.GetFFT_MaxUnitCell_A() > indexing_bound_A) { - IndexerThreadPool wide_pool(wide_settings, IndexerConstruction::OnFirstUse); - constrained = pick_best(wide_pool, *indexer); - } else { - constrained = pick_best(*indexer_pool, *indexer); - } - experiment_.SetUnitCell(std::nullopt); // leave the space-group / cell determination de-novo - if (constrained.result.has_value() && constrained.score > best.score) - best = constrained; - } - } + // The long-axis rescue (above), at the centre the run now stands on. + if (!long_axis_asked) + if (auto rescued = long_axis_rescue(best)) + best = *rescued; release_spot_engines(); // no more frames to read; give the cards their memory back if (!cancelled_) {