diff --git a/docs/RUGNUX.md b/docs/RUGNUX.md index 3df79800..a1e3651d 100644 --- a/docs/RUGNUX.md +++ b/docs/RUGNUX.md @@ -34,9 +34,11 @@ scale with the thread count (`-N`). ## Input and output -**Input** is a single Jungfraujoch HDF5 master file (NXmx-based). If the dataset already contains -stored spot lists, two-pass rotation indexing can reuse them instead of re-running spot finding on -the first pass. +**Input** is a single Jungfraujoch HDF5 master file (NXmx-based). Spots are always found by `rugnux` +itself, including for the two-pass rotation first pass — the spot lists a dataset may already carry +were found online, at the acquisition's threshold and with its ice-band spots already discarded, so +reusing them would hide the spot-finding settings from the lattice search and stop the run measuring +its own ice. **Output** (controlled by `-o, --output-prefix`, default `output`): @@ -240,7 +242,6 @@ rotation explicitly and pick the pass or lattice. | `-r, --refine ` | Geometry refinement: `none` \| `orientation` \| `beam_and_lattice` (default) \| `flex` (try all three per image, keep whichever indexes the most spots; alias `multi`) | | `-R, --two-pass-rotation[=num]` | Two-pass offline rotation indexing (default for goniometer data; optional first-pass image count, default 100) | | `--single-pass-rotation[=num]` | Online-like single-pass rotation indexing (optional min angular range, deg) | -| `--redo-rotation-spots` | Redo spot finding for the two-pass rotation first pass | | `--force-rotation-lattice ` | Force rotation lattice (9 floats, Å), skipping the first pass | | `--rotation-no-postrefine` | Rotation: disable the default-on two-pass geometry post-refine (see the rotation section) | | `--refine-geometry[=N\|off]` | Stills: extra first pass that bundle-adjusts the shared beam/distance/cell from N strongly-indexed frames (default 200) then re-indexes; default ON for stills with a reference cell (`-C` / `-z`), `=off` disables | diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 0b979636..502a8b75 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -64,8 +64,7 @@ namespace { IndexAndRefine::IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool *indexer, bool retain_outcomes, bool real_time) - : index_ice_rings(x.GetIndexingSettings().GetIndexIceRings()), - retain_outcomes_(retain_outcomes), + : retain_outcomes_(retain_outcomes), real_time(real_time), experiment(x), geom_(x.GetDiffractionGeometry()), @@ -158,6 +157,9 @@ IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetry(Data IndexerResult indexer_result; bool any_executed = false; float best_frac = -1.0f; + // Read per call rather than cached at construction: the run measures for itself whether the + // crystal has ice, and that verdict lands after this object is built. + const bool index_ice_rings = experiment.GetIndexingSettings().GetIndexIceRings(); for (size_t seed_cap : {size_t{30}, size_t{80}, std::numeric_limits::max()}) { std::vector recip; recip.reserve(std::min(seed_cap, msg.spots.size())); diff --git a/image_analysis/IndexAndRefine.h b/image_analysis/IndexAndRefine.h index 7c9268de..d79c2b0a 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -29,7 +29,6 @@ using BraggIntegrateFn = std::function( const std::vector &predicted, size_t npredicted, int64_t image_number)>; class IndexAndRefine { - const bool index_ice_rings; // When false, the current image's result is still returned via the outgoing message, but the // whole-run integration_outcome vector is not retained (viewer live/interactive use, which never // scales the accumulated run). rugnux/receiver keep it true so ScaleAllImages/merge have the data. diff --git a/image_analysis/spot_finding/SpotUtils.cpp b/image_analysis/spot_finding/SpotUtils.cpp index 304653cf..641ebd26 100644 --- a/image_analysis/spot_finding/SpotUtils.cpp +++ b/image_analysis/spot_finding/SpotUtils.cpp @@ -70,14 +70,17 @@ void MarkIceRings(std::vector &spots, float tolerance_q_recipA) { } } -void FilterSpotsByCount(std::vector &input, int64_t count) { +void FilterSpotsByCount(std::vector &input, int64_t count, bool deprioritise_ice) { size_t output_size = std::min(input.size(), count); std::ranges::partial_sort(input, input.begin() + output_size, std::ranges::less{}, // comparator on the projected key - [](const SpotToSave &s) { - // projection: non-ice first (false < true), then strongest intensity first. - return std::tuple{s.ice_ring, -s.intensity}; + [deprioritise_ice](const SpotToSave &s) { + // projection: non-ice first (false < true), then strongest intensity + // first. Where the run has no measurable ice the flag marks ordinary + // reflections that happen to lie in the fixed bands, so ordering on it + // would discard a fifth of the strongest spots for nothing. + return std::tuple{deprioritise_ice && s.ice_ring, -s.intensity}; }); input.resize(output_size); } @@ -196,7 +199,10 @@ void SpotAnalyze(const DiffractionExperiment &experiment, output.resolution_estimate = GetResolution(spots_out); - FilterSpotsByCount(spots_out, experiment.GetMaxSpotCount()); + // One decision drives both: if indexing is to use the ice-band spots, the spot budget must not + // throw them away before it gets the chance. + FilterSpotsByCount(spots_out, experiment.GetMaxSpotCount(), + !experiment.GetIndexingSettings().GetIndexIceRings()); output.spots = spots_out; } diff --git a/image_analysis/spot_finding/SpotUtils.h b/image_analysis/spot_finding/SpotUtils.h index 23ca156d..a6b79e42 100644 --- a/image_analysis/spot_finding/SpotUtils.h +++ b/image_analysis/spot_finding/SpotUtils.h @@ -20,7 +20,10 @@ float CountIceRingControlSpots(const std::vector &spots, float half_ void MarkIceRings(std::vector &spots, float tolerance_q_recipA); -void FilterSpotsByCount(std::vector &input, int64_t count); +// Keep the strongest `count` spots. With deprioritise_ice, spots on the hexagonal ice bands are ordered +// last and so are the first to go; pass false where the run has no measurable ice, in which case the +// flag marks ordinary reflections and ordering on it would discard good data. +void FilterSpotsByCount(std::vector &input, int64_t count, bool deprioritise_ice); void FilterSpuriousHighResolutionSpots(std::vector &spots, float threshold); // Ignore high res. spots if there is a gap in (1/d) between two spots of dist_threshold (default: 0.25 A^-1) diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 94d3abf9..00efab99 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -563,6 +563,11 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b const auto start_time = std::chrono::steady_clock::now(); // First pass of two-pass rotation indexing (full analysis only). + // Ice measured on the first-pass sample (see where they are filled): pooled here so the verdict + // below, which runs after the first pass and before the per-image loop, can read them. + std::vector sample_ice_scores; + double sample_ice_ring_spots = 0.0, sample_ice_control_spots = 0.0; + if (full && force_rotation_result_.has_value()) { // Supercell-collapse fallback: force pass-1's WHOLE indexing result (lattice + refined orientation / // search metadata / axis), not just its lattice - a lattice-only force loses that metadata and @@ -583,36 +588,17 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b return std::nullopt; }; - // Reuse stored spots only if the file actually has them: a plain DECTRIS dataset carries no - // spot-finding results, so "reusing" would leave every frame with zero spots and fail - // indexing - find them instead. (--redo-rotation-spots already forces finding.) - const bool reuse_spots = config_.reuse_rotation_spots && reader_.HasStoredSpots(); - if (config_.reuse_rotation_spots && !reuse_spots) - logger.Info("Dataset has no stored spots - running spot finding for first-pass rotation indexing"); - // Say so when the stored spots are used: the lattice is then determined - and cross-validated - - // from the spots the acquisition wrote, not from anything the spot-finding settings produce, so a - // reader comparing settings needs to know this pass did not see them. - if (reuse_spots) - logger.Info("Using the spots stored in the file for first-pass rotation indexing " - "(spot-finding settings apply to the per-image pass only; --redo-rotation-spots " - "re-finds them here too)"); - // Spots for a first-pass image ordinal, cached (a frame is reused across schemes and the // validation set; re-finding is expensive on the --redo-rotation-spots path). Accessed only // from single-threaded sections (the scheme feed and the validation loop), so no locking. std::map> spot_cache; // One analysis engine for the whole first pass. Constructing it allocates a CUDA stream and a // full set of GPU buffers (preprocessing, spot finding, azimuthal integration, Bragg - // integration) - far more work than analysing a frame - so building it per image made - // --redo-rotation-spots pay for hundreds of them, serially. Only needed when spots have to be - // found; this section is single-threaded, so one engine is enough. - std::unique_ptr analysis; - std::unique_ptr profile; - if (!reuse_spots) { - analysis = std::make_unique(experiment_, mapping, pixel_mask_, *indexer, - /*enable_fused_adaptive_gpu=*/true); - profile = std::make_unique(mapping); - } + // integration) - far more work than analysing a frame - so building it per image would pay for + // hundreds of them, serially. This section is single-threaded, so one engine is enough. + auto analysis = std::make_unique(experiment_, mapping, pixel_mask_, *indexer, + /*enable_fused_adaptive_gpu=*/true); + auto profile = std::make_unique(mapping); auto get_spots = [&](int ordinal) -> const std::vector & { auto it = spot_cache.find(ordinal); if (it != spot_cache.end()) @@ -620,9 +606,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b const int image_idx = start_image + ordinal * config_.stride; std::vector spots; try { - if (reuse_spots) { - spots = reader_.ReadSpots(image_idx); - } else if (auto img = reader_.GetRawImage(image_idx)) { + if (auto img = reader_.GetRawImage(image_idx)) { DataMessage m{}; m.number = ordinal; m.original_number = image_idx; @@ -633,6 +617,12 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b if (dataset->efficiency.size() > image_idx) m.image_collection_efficiency = dataset->efficiency[image_idx]; analysis->Analyze(m, *profile, first_pass); + if (m.ice_ring_score) + sample_ice_scores.push_back(*m.ice_ring_score); + if (m.spot_count_ice_rings) + sample_ice_ring_spots += static_cast(*m.spot_count_ice_rings); + if (m.spot_count_ice_control) + sample_ice_control_spots += static_cast(*m.spot_count_ice_control); spots = std::move(m.spots); } } catch (const std::exception &e) { @@ -889,6 +879,39 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b } } + // Ice-presence verdict, from the first-pass sample and before the pass that acts on it. The eleven + // hexagonal bands are fixed geometry and hold 16-26 % of the unique reflections whether or not the + // crystal has ice, so setting their spots aside on a clean crystal discards a fifth of them for + // nothing - and de-prioritising them in the --max-spots budget throws away the strongest first. + // Where the sample sees no ice, index on all the spots instead. Only the SAMPLE can decide this: + // the spot channel is a ratio pooled over images, so it is meaningless on any single frame. + if (full && !sample_ice_scores.empty() && experiment_.IsDetectIceRings() + && !experiment_.GetIndexingSettings().GetIndexIceRings()) { + const float min_score = experiment_.GetScalingSettings().GetIceMinScore(); + const float min_ratio = experiment_.GetScalingSettings().GetIceMinSpotRatio(); + const float score = std::accumulate(sample_ice_scores.begin(), sample_ice_scores.end(), 0.0f) + / static_cast(sample_ice_scores.size()); + const bool smooth = score >= min_score; + const bool spotty = min_ratio > 0.0f && sample_ice_control_spots > 0.0 + && sample_ice_ring_spots / sample_ice_control_spots >= min_ratio; + if (!smooth && !spotty) { + IndexingSettings is = experiment_.GetIndexingSettings(); + is.IndexIceRings(true); + experiment_.ImportIndexingSettings(is); + logger.Info("First-pass sample of {} images sees no ice (score {:.2f} < {:.2f}, spot ratio " + "{:.2f} < {:.2f}): indexing on the ice-band spots too", + sample_ice_scores.size(), score, min_score, + sample_ice_control_spots > 0.0 + ? sample_ice_ring_spots / sample_ice_control_spots : NAN, min_ratio); + } else { + logger.Info("First-pass sample of {} images sees ice (score {:.2f}, spot ratio {:.2f}): " + "ice-band spots set aside for indexing", + sample_ice_scores.size(), score, + sample_ice_control_spots > 0.0 + ? sample_ice_ring_spots / sample_ice_control_spots : NAN); + } + } + // Main per-image loop, spread over N worker threads pulling from a shared counter. HDF5 reads // are serialized by the global hdf5_mutex; the analysis runs in parallel. std::atomic next_ordinal = 0; diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 7d524d7a..7da26936 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -59,7 +59,6 @@ struct ProcessConfig { // Rotation indexing (FullAnalysis) bool rotation_indexing = false; bool two_pass_rotation = true; - bool reuse_rotation_spots = true; int rotation_indexing_image_count = 100; std::optional forced_rotation_lattice; @@ -200,6 +199,7 @@ class Rugnux { // the search chose, which is the circular claim that reasoning exists to flag. bool prepass_promoted_point_group_ = false; + // Stills global geometry-refinement first pass (config_.refine_geometry): index a spread sample of // frames, bundle-adjust the shared beam/distance/cell from the strongest ones, and apply the result // to experiment_ so the main pass re-indexes with it. No-op (leaves experiment_ unchanged) if it diff --git a/rugnux/RugnuxCommandLine.cpp b/rugnux/RugnuxCommandLine.cpp index 91d5769d..16dec5ad 100644 --- a/rugnux/RugnuxCommandLine.cpp +++ b/rugnux/RugnuxCommandLine.cpp @@ -141,8 +141,6 @@ std::string RugnuxCommandLine(const ProcessConfig &config, args.push_back("-R" + std::to_string(config.rotation_indexing_image_count)); else args.emplace_back("--single-pass-rotation"); - if (!config.reuse_rotation_spots) - args.emplace_back("--redo-rotation-spots"); } else if (experiment.GetGoniometer().has_value()) { // rotation dataset processed as stills -> the user overrode the default with --force-still args.emplace_back("--force-still"); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 93cc6914..77fbc64c 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -86,7 +86,6 @@ void print_usage() { std::cout << " --force-still Process a rotation (goniometer) dataset as independent stills (still indexing + per-image ScaleOnTheFly) instead of rotation" << std::endl; std::cout << " -R, --two-pass-rotation[=num] Two-pass offline rotation indexing (default for goniometer data; optional first-pass image count, default: 100)" << std::endl; std::cout << " --single-pass-rotation[=num] Use online-like single-pass rotation indexing (optional: min angular range deg)" << std::endl; - std::cout << " --redo-rotation-spots Redo spot finding for two-pass rotation indexing" << std::endl; std::cout << " --force-rotation-lattice Force rotation indexer with external lattice (in Angstrom) : \"a0x,a0y,a0z,a1x,a1y,a1z,a2x,a2y,a2z\" (9 floats, skips first pass)" << std::endl; std::cout << " --rotation-no-postrefine Disable the (default-on) two-pass rotation post-refine (post-refine detector distance/beam + cell/axis, then re-integrate; the refined pass is the canonical _* output, the header-geometry pass is kept as _01_*)" << std::endl; std::cout << " -X, --indexing-algorithm Indexing algorithm (FFBIDX|FFT|FFTW|Auto|None)" << std::endl; @@ -173,7 +172,6 @@ enum { OPT_RESOLUTION_CC_TARGET, OPT_RESOLUTION_SHELLS, OPT_SINGLE_PASS_ROTATION, - OPT_REDO_ROTATION_SPOTS, OPT_FORCE_ROTATION_LATTICE, OPT_ROTATION_NO_POSTREFINE, OPT_BACKGROUND_CLIP, @@ -265,7 +263,6 @@ static option long_options[] = { {"rot1", required_argument, nullptr, OPT_ROT1}, {"rot2", required_argument, nullptr, OPT_ROT2}, {"polarization", required_argument, nullptr, OPT_POLARIZATION}, - {"redo-rotation-spots", no_argument, nullptr, OPT_REDO_ROTATION_SPOTS}, {"force-rotation-lattice", required_argument, nullptr, OPT_FORCE_ROTATION_LATTICE}, {"rotation-no-postrefine", no_argument, nullptr, OPT_ROTATION_NO_POSTREFINE}, {"refine-geometry", optional_argument, nullptr, OPT_REFINE_GEOMETRY}, @@ -524,11 +521,9 @@ static int RunRugnux(int argc, char **argv) { bool rotation_indexing = false; bool force_still = false; // --force-still: process a rotation dataset as stills (indexing + scaling) bool two_pass_rotation = true; - bool reuse_rotation_spots = true; // Set by any spot-finding option. The two-pass rotation first pass reuses the spots stored in the // file when it has them, so those options would otherwise not reach the pass that determines the // lattice - the setting would appear to do nothing at all on rotation data. - bool spot_finding_options_given = false; bool rotation_postrefine_geometry = true; // default on; --rotation-no-postrefine disables it int rotation_indexing_image_count = 100; std::optional rotation_indexing_range; @@ -647,9 +642,6 @@ static int RunRugnux(int argc, char **argv) { if (optarg) rotation_indexing_range = atof(optarg); break; - case OPT_REDO_ROTATION_SPOTS: - reuse_rotation_spots = false; - break; case OPT_ROTATION_NO_POSTREFINE: rotation_postrefine_geometry = false; break; @@ -777,45 +769,37 @@ static int RunRugnux(int argc, char **argv) { break; } case OPT_SPOT_SIGMA: - spot_finding_options_given = true; sigma_spot_finding = parse_number_arg(optarg, "--spot-sigma", logger, 1.0f); logger.Info("Noise threshold level for spot finding set to {:.2f} sigma", sigma_spot_finding); break; case OPT_SPOT_THRESHOLD: - spot_finding_options_given = true; photon_count_threshold_spot_finding = parse_number_arg(optarg, "--spot-threshold", logger, 0); logger.Info("Photon-count threshold level for spot finding set to {:d}", photon_count_threshold_spot_finding); break; case OPT_MIN_PIX_PER_SPOT: - spot_finding_options_given = true; // Giving an explicit min-pix opts out of the per-image adaptive selection. min_pix_per_spot = parse_number_arg(optarg, "--min-pix-per-spot", logger, 1); logger.Info("Minimum pixels per spot fixed at {:d} (adaptive per-image min-pix off)", *min_pix_per_spot); break; case OPT_ADAPTIVE_SPOTS: - spot_finding_options_given = true; adaptive_spots = true; logger.Info("Adaptive (self-calibrating) spot detection enabled"); break; case OPT_NO_ADAPTIVE_SPOTS: - spot_finding_options_given = true; adaptive_spots = false; logger.Info("Adaptive spot detection off: using the fixed --spot-threshold / --spot-sigma finder"); break; case OPT_SPOT_FALSE_PIXELS: - spot_finding_options_given = true; false_pixels_per_frame = parse_number_arg(optarg, "--spot-false-pixels", logger, 1.0f); adaptive_spots = true; logger.Info("Adaptive spot detection: expected false pixels/frame set to {:.0f}", false_pixels_per_frame); break; case OPT_SPOT_LOW_RESOLUTION: - spot_finding_options_given = true; d_max_spot_finding = parse_number_arg(optarg, "--spot-low-resolution", logger, 0.0f); logger.Info("Low resolution limit for spot finding set to {:.1f} A", d_max_spot_finding); break; case OPT_SPOT_RESOLUTION: { - spot_finding_options_given = true; // 0 has always meant "no limit" for this setting; keep that, but express it as the unset // optional the rest of the code understands. Passing the 0 through instead reached // ResolutionShells (via the spot plot), which rejects a zero d_min and threw away every image. @@ -830,7 +814,6 @@ static int RunRugnux(int argc, char **argv) { break; } case OPT_MAX_SPOTS: - spot_finding_options_given = true; max_spot_count_override = parse_number_arg(optarg, "--max-spots", logger, 1); break; case OPT_AZINT_ONLY: @@ -849,8 +832,6 @@ static int RunRugnux(int argc, char **argv) { scale_fulls_arg = false; break; case OPT_DETECT_ICE_RINGS: - // Deliberately NOT a spot-finding option for the purpose of re-finding the first-pass - // rotation spots - see where reuse_rotation_spots is decided below. if (optarg == nullptr || strcmp(optarg, "on") == 0) detect_ice_rings = true; else if (strcmp(optarg, "off") == 0) @@ -1592,27 +1573,6 @@ static int RunRugnux(int argc, char **argv) { else if (!dataset->file_detect_ice_rings.has_value()) experiment.DetectIceRings(rotation_indexing); - // First-pass rotation indexing reuses the spots stored in the file, and does NOT re-mark them - - // their ice flags are the ones the acquisition wrote. So --detect-ice-rings can only invalidate - // them when it asks for something the file did not do; matching the file is a no-op for that pass. - // Re-finding whenever the flag is merely present would swap the acquisition's spots for this - // program's own, which moves the first-pass lattice by itself: measured over the rotation battery, - // passing the semantically null --detect-ice-rings=on to files that already carry it changed the - // merged data on every crystal and lost one to indexing failure. That also made the flag impossible - // to test, since both arms of an A/B moved for a reason unrelated to ice. - // A file with no key at all is treated as "did not mark", which is what its stored spots show: - // such a dataset carries no per-spot ice flags to reuse. That is a statement about the SPOTS, not - // about the setting - which is why it stays value_or(false) even though a keyless rotation file - // now defaults to detecting ice. - if (detect_ice_rings.has_value() - && detect_ice_rings.value() != dataset->file_detect_ice_rings.value_or(false)) { - if (reuse_rotation_spots) - logger.Info("--detect-ice-rings={} differs from the spots stored in the file: re-finding " - "them for first-pass rotation indexing too", - detect_ice_rings.value() ? "on" : "off"); - reuse_rotation_spots = false; - } - // Scale-fulls refits the per-frame scale on the rotation combined fulls; on by default for rotation // data (where it lifts ISa substantially) and off for stills. --no-scale-fulls overrides. const bool scale_fulls = scale_fulls_arg.value_or(rotation_indexing); @@ -1787,15 +1747,6 @@ static int RunRugnux(int argc, char **argv) { config.spot_finding = spot_settings; config.rotation_indexing = rotation_indexing; config.two_pass_rotation = two_pass_rotation; - // Asking for particular spots means asking the lattice search to use them. Without this the stored - // spots win on any file that has them, so every spot-finding option is a no-op for rotation indexing - // and the settings appear to have no effect whatsoever on the cell that comes out. - if (spot_finding_options_given && reuse_rotation_spots) { - reuse_rotation_spots = false; - logger.Info("Spot-finding settings given: re-finding spots for first-pass rotation indexing too, " - "so the lattice search sees them"); - } - config.reuse_rotation_spots = reuse_rotation_spots; config.rotation_postrefine_geometry = rotation_postrefine_geometry; config.rotation_indexing_image_count = rotation_indexing_image_count; config.forced_rotation_lattice = forced_rotation_lattice; diff --git a/tests/RugnuxLargeTest.cpp b/tests/RugnuxLargeTest.cpp index 1b695641..d8b99587 100644 --- a/tests/RugnuxLargeTest.cpp +++ b/tests/RugnuxLargeTest.cpp @@ -65,7 +65,6 @@ TEST_CASE("Rugnux_Rotation", "[large]") { config.spot_finding.indexing = true; config.rotation_indexing = true; config.two_pass_rotation = true; - config.reuse_rotation_spots = false; // redo spot finding (raw dataset may carry no spots) Rugnux process(reader, experiment, *dataset->pixel_mask, config); ProcessResult result;