diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 8efde9269..68b55bd77 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -869,7 +869,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b std::future saturation; if (full && !experiment_.IsRotationIndexing() && !config_.live) saturation = std::async(std::launch::async, [&] { return SampleSaturation(start_image, images_to_process); }); - if (full && config_.refine_geometry.has_value()) + if (full && config_.refine_geometry.has_value() && !config_.detector_distance_unknown) RefineStillsGeometry(start_image, end_image, images_to_process, observer); // The geometry pre-pass integrates only the first part of the sweep (pass 2 reads it all), and on a diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 63ea28467..bdb296d91 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -177,9 +177,12 @@ struct ProcessConfig { // two answers - the one comparison an on-failure trigger cannot make, because a centre wrong // along the spindle does not fail. bool beam_center_check = true; - // The detector distance was given on the command line: stills do not search it (DecideStillsGeometry); - // the stills geometry refinement starts from it. + // The detector distance was given on the command line: on stills it is held - not searched + // (DecideStillsGeometry), not refined (RefineStillsGeometry). bool detector_distance_given = false; + // Neither the input nor the command line states the distance (an XTC): the input's value is a + // placeholder, which the stills search does not hold against the others and nothing refines from. + bool detector_distance_unknown = false; // Whether a saturation level the frames show below the file's limit replaces it (DetectSaturation): // for every input but a Jungfraujoch-written file, whose stated limit is trusted. bool detect_saturation = false; diff --git a/rugnux/RugnuxStills.cpp b/rugnux/RugnuxStills.cpp index ac0e8e139..af3849fbe 100644 --- a/rugnux/RugnuxStills.cpp +++ b/rugnux/RugnuxStills.cpp @@ -100,6 +100,55 @@ void AssignModules(std::vector &frames, const DetectorModules & s.module = static_cast(modules.ModuleOfPlanePoint(s.x, s.y)); } +// A frame the geometry refinement can use: its spots on lattice 0 with their index, the others without. +GeomRefineFrame ToRefineFrame(const DataMessage &msg, int ordinal) { + GeomRefineFrame f; + f.lattice = *msg.indexing_lattice; + f.ordinal = ordinal; + for (const auto &s : msg.spots) { + if (s.indexed && s.lattice == 0) + f.spots.push_back(GeomRefineSpot{s.x, s.y, + static_cast(s.h), static_cast(s.k), static_cast(s.l), s.panel}); + else + f.unindexed.push_back(GeomRefineSpot{s.x, s.y, 0, 0, 0, s.panel}); + } + return f; +} + +// Spots piled up at one radius - an ice or other powder ring - are not this crystal's lattice. They are +// found by their radius in millimetres about (beam_x, beam_y), not by a d-spacing: a d computed at a +// distance that may be wrong misses the rings, and then they are most of the spots a frame has. +// A 0.5 mm annulus is a ring where it holds at least 10 spots and over five times the median of the +// annuli within 10 mm of it; the annuli either side go with it. Returns how many spots were dropped. +size_t DropPowderRingSpots(std::vector &frames, float beam_x, float beam_y, float pixel_mm) { + constexpr float BIN_MM = 0.5f; + constexpr int WINDOW_BINS = 20, MIN_SPOTS = 10, OVER_MEDIAN = 5; + const auto bin_of = [&](const SpotToSave &s) { + return static_cast(std::hypot(s.x - beam_x, s.y - beam_y) * pixel_mm / BIN_MM); + }; + std::vector count; + for (const auto &f : frames) + for (const auto &s : f.spots) { + const size_t b = bin_of(s); + if (b >= count.size()) + count.resize(b + 1, 0); + count[b]++; + } + std::vector ring(count.size(), 0); + for (size_t i = 0; i < count.size(); i++) { + std::vector w(count.begin() + static_cast(i > WINDOW_BINS ? i - WINDOW_BINS : 0), + count.begin() + static_cast(std::min(count.size(), i + WINDOW_BINS + 1))); + std::nth_element(w.begin(), w.begin() + w.size() / 2, w.end()); + if (count[i] >= MIN_SPOTS && count[i] > OVER_MEDIAN * std::max(w[w.size() / 2], 1)) + for (size_t j = (i > 0 ? i - 1 : 0); j <= std::min(count.size() - 1, i + 1); j++) + ring[j] = 1; + } + size_t dropped = 0; + for (auto &f : frames) + dropped += std::erase_if(f.spots, [&](const SpotToSave &s) { return ring[bin_of(s)] != 0; }); + return dropped; +} + } // namespace void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_process, @@ -214,16 +263,7 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_ if (!msg.indexing_result.value_or(false) || !msg.indexing_lattice.has_value()) continue; - GeomRefineFrame f; - f.lattice = *msg.indexing_lattice; - f.ordinal = ordinal; - for (const auto &s : msg.spots) { - if (s.indexed && s.lattice == 0) - f.spots.push_back(GeomRefineSpot{s.x, s.y, - static_cast(s.h), static_cast(s.k), static_cast(s.l), s.panel}); - else - f.unindexed.push_back(GeomRefineSpot{s.x, s.y, 0, 0, 0, s.panel}); - } + GeomRefineFrame f = ToRefineFrame(msg, ordinal); if (f.spots.size() >= 6) { if (static_cast(f.spots.size()) >= MIN_STRONG_SPOTS) stripe_strong++; @@ -295,6 +335,64 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_ return; } + // The bundle's spots were indexed at the geometry the fit started from, and the high-angle ones only + // where that geometry put them close enough to their reflections: a fit that has moved the distance far + // carries the start in which spots it rests on (two starts a percent apart ended 0.2 % apart). So where it moved the distance by more than REFIT_ABOVE, the bundle's frames are indexed + // again at the fitted geometry and cell and fitted again, until the distance stays put. Not where + // modules are being placed: the placement below is fitted on spots indexed at the start, as before. + constexpr double REFIT_ABOVE = 0.003, REFIT_CONVERGED = 0.0005; + constexpr int MAX_REFITS = 3; + double moved_by = std::abs(r.distance_mm / nominal.GetDetectorDistance_mm() - 1.0); + for (int refit = 0; refit < MAX_REFITS && !placing_modules && moved_by > REFIT_CONVERGED + && (refit > 0 || moved_by > REFIT_ABOVE) && !cancelled_; refit++) { + TimingMark refit_mark(logger, "geometry refit"); + DiffractionExperiment e = experiment_; + e.BeamX_pxl(r.beam_x_px).BeamY_pxl(r.beam_y_px).DetectorDistance_mm(r.distance_mm); + e.SetUnitCell(r.cell); + AzimuthalIntegrationMapping refit_mapping(e, pixel_mask_); + IndexerThreadPool refit_pool(e.GetIndexingSettings(), IndexerConstruction::OnFirstUse, gpu_work_beside_); + IndexAndRefine refit_indexer(e, &refit_pool, /*retain_outcomes=*/false); + std::vector again(strong.size()); + std::atomic next = 0; + auto refit_worker = [&]() { + pin_gpu(); + JFJochReaderRawImage img; + MXAnalysisWithoutFPGA analysis(e, refit_mapping, pixel_mask_, refit_indexer, + /*enable_fused_adaptive_gpu=*/true); + AzimuthalIntegrationProfile profile(refit_mapping); + for (size_t i = next.fetch_add(1); i < strong.size() && !cancelled_; i = next.fetch_add(1)) { + const int image_idx = start_image + strong[i].ordinal * config_.stride; + if (!reader_.ReadRawImage(image_idx, img)) + continue; + DataMessage msg{}; + msg.image = img.image; + msg.number = strong[i].ordinal; + msg.original_number = image_idx; + if (dataset->efficiency.size() > image_idx) + msg.image_collection_efficiency = dataset->efficiency[image_idx]; + analysis.Analyze(msg, profile, pass_settings); + if (msg.indexing_result.value_or(false) && msg.indexing_lattice.has_value()) + again[i] = ToRefineFrame(msg, strong[i].ordinal); + } + }; + std::vector> refit_futures; + for (int i = 0; i < ImageWorkers(); ++i) + refit_futures.push_back(std::async(std::launch::async, refit_worker)); + for (auto &fut : refit_futures) + fut.get(); + refit_indexer.RethrowIndexerError(); + std::erase_if(again, [&](const GeomRefineFrame &f) { + return static_cast(f.spots.size()) < MIN_STRONG_SPOTS; + }); + const GeometryRefinerResult r2 = RefineGlobalGeometry(e.GetDiffractionGeometry(), r.cell, again, gr_settings); + if (!r2.ok) + break; + moved_by = std::abs(r2.distance_mm / r.distance_mm - 1.0); + logger.Info("Geometry refinement: indexed again at the fitted geometry ({} frames), distance {:.4f} -> {:.4f} mm", + again.size(), r.distance_mm, r2.distance_mm); + r = r2; + } + // A detector of several modules: each placed as a rigid block, where that predicts the spots better on // frames the placement was not fitted to. The strong frames are dealt alternately into two halves; each // half is fitted with the modules held and with them moved, and each fit is scored on the other half @@ -431,6 +529,11 @@ std::vector Rugnux::FindStillsSampleSpots(const std::vector &o settings.enable = true; settings.indexing = false; settings.quick_integration = false; + // No resolution limits: they are d-spacings, evaluated at the distance the run has now, and this sample + // is what decides whether that distance is right. At a file's 839 mm for a detector at 106 mm the 50 A + // limit alone cut every spot within 20 mm of the beam. + settings.low_resolution_limit = std::nullopt; + settings.high_resolution_limit = std::nullopt; std::vector frames(ordinals.size()); std::atomic next = 0; @@ -520,27 +623,45 @@ bool Rugnux::DecideStillsGeometry(int start_image, int images_to_process) { TimingMark mark(logger, "stills geometry trial"); constexpr int SAMPLE_FRAMES = 96; - // The distance grid is scanned on the frames with the most spots only: a frame with few has little - // at high angle, where the distance shows. - constexpr int GRID_FRAMES = 32; + // The distance is told apart far from the beam, so its searches run on the frames with the most spots: + // the scan below on SCAN_FRAMES of them (it only has to find the peak), the steps about each peak on + // GRID_FRAMES. + constexpr int SCAN_FRAMES = 8, GRID_FRAMES = 32; // Below this fraction of the best hypothesis' count the file's value is outside the peak. constexpr double FILE_KEPT_ABOVE = 0.9; - // The distance grid: the file's value times STEP^k, over [MIN_FACTOR, MAX_FACTOR]. Wide enough for an - // encoder reading or a nominal camera length to be wrong, fine enough not to step over the peak. - constexpr float STEP = 1.025f, MIN_FACTOR = 0.6f, MAX_FACTOR = 1.6f; - - const std::vector sample = select_equally_spaced_image_ordinals( - images_to_process, std::min(images_to_process, SAMPLE_FRAMES)); - std::vector frames = FindStillsSampleSpots(sample, start_image); - std::stable_sort(frames.begin(), frames.end(), [](const DataMessage &a, const DataMessage &b) { - return a.spots.size() > b.spots.size(); - }); - const std::vector strongest(frames.begin(), - frames.begin() + std::min(frames.size(), GRID_FRAMES)); - IndexerThreadPool pool(experiment_.GetIndexingSettings(), IndexerConstruction::OnFirstUse, gpu_work_beside_); + // The distance scan: STEP^k mm over [MIN_DISTANCE_MM, MAX_DISTANCE_MM]. The range is in millimetres, + // not in multiples of the file's value, so where the scan looks does not depend on the value it checks + // (a Cheetah file states 839 mm for a detector at 106 mm). The step is the peak's width: a distance + // 1.2 % off the peak still puts a third of its spots on the lattice, 3.6 % off almost none. + constexpr float STEP = 1.025f, MIN_DISTANCE_MM = 30.0f, MAX_DISTANCE_MM = 3000.0f; + // How many of the scan's local maxima are looked at more closely. + constexpr size_t SCAN_PEAKS = 3; + // A sample found at a distance this far from the one decided is found again there (see below). + constexpr float RESAMPLE_ABOVE = 0.05f; const float file_x = experiment_.GetBeamX_pxl(), file_y = experiment_.GetBeamY_pxl(); const float file_distance = experiment_.GetDetectorDistance_mm(); + const std::vector sample = select_equally_spaced_image_ordinals( + images_to_process, std::min(images_to_process, SAMPLE_FRAMES)); + std::vector frames, strongest, scan_frames; + // The sample's spots, found at the run's present geometry; powder rings off about the background's + // centre where there is one (it does not depend on the distance), strongest frames first. + const auto find_sample = [&]() { + frames = FindStillsSampleSpots(sample, start_image); + const float bx = background_center_ ? background_center_->beam_x_pxl : experiment_.GetBeamX_pxl(); + const float by = background_center_ ? background_center_->beam_y_pxl : experiment_.GetBeamY_pxl(); + const size_t dropped = DropPowderRingSpots(frames, bx, by, experiment_.GetPixelSize_mm()); + if (dropped > 0) + logger.Info("Stills geometry trial: {} sampled spots on powder rings left out", dropped); + std::stable_sort(frames.begin(), frames.end(), [](const DataMessage &a, const DataMessage &b) { + return a.spots.size() > b.spots.size(); + }); + strongest.assign(frames.begin(), frames.begin() + std::min(frames.size(), GRID_FRAMES)); + scan_frames.assign(frames.begin(), frames.begin() + std::min(frames.size(), SCAN_FRAMES)); + }; + find_sample(); + IndexerThreadPool pool(experiment_.GetIndexingSettings(), IndexerConstruction::OnFirstUse, gpu_work_beside_); + // A hypothesis is indexed once: the searches below meet the same one more than once - the decided centre // is the file's whenever the file's is kept - and it gives the same count every time. std::map *>, int64_t> tried; @@ -590,45 +711,68 @@ bool Rugnux::DecideStillsGeometry(int start_image, int images_to_process) { score_here = moved ? best : file_score; }; decide_centre(file_distance); + // A file that states no distance (an XTC) has a placeholder there, not a hypothesis. + if (config_.detector_distance_unknown) + score_here = 0; - // The distance, where a reference cell is given: a grid at the centre just decided, on the strongest - // frames, then a finer one about its best on the whole sample, which also decides against the file's - // value. With the cell fixed each frame is matched rather than searched for, so the grid is cheap - - // and only then is the distance a hypothesis of its own: without a cell a wrong distance indexes as a - // cell of the wrong size, which only the curvature far from the beam tells apart. The stills geometry - // refinement takes it from there (RunPipeline). - // A distance the user gave on the command line is not searched: a sample of mixed crystal forms can - // put more spots on one given cell at a distance that is not the detector's. Like any other, it is - // where the stills geometry refinement starts, not a value it holds. + // The distance, where a reference cell is given. With the cell fixed each frame is matched rather than + // searched for, so a trial is cheap - and only then is the distance a hypothesis of its own: without a + // cell a wrong distance indexes as a cell of the wrong size, which only the curvature far from the beam + // tells apart. The stills geometry refinement takes it from there (RunPipeline). + // + // Like the beam centre's capture, the search is global first and local after: a scan over the whole + // range at one centre that does not depend on the file (the background's, else the capture's), then + // each of its strongest peaks a step either side at every centre hypothesis, then a finer step about the + // best on the whole sample, which also decides against the file's value. A distance the user gave + // (--detector-distance) is not searched: it is the user's decision, e.g. a calibrated distance where + // the reference cell is the less certain of the two. float distance = file_distance; if (experiment_.GetUnitCell().has_value() && !config_.detector_distance_given) { - // Every centre hypothesis gets its own grid: where the file's distance is wrong, the centre - // decided at it had little indexing to go on, and a wrong centre in turn leaves the grid flat. + std::vector grid; + for (float d = MIN_DISTANCE_MM; d <= MAX_DISTANCE_MM; d *= STEP) + grid.push_back(d); + const std::pair scan_centre = !candidates.empty() ? candidates.front() : std::pair{x, y}; + std::vector scan(grid.size()); + for (size_t i = 0; i < grid.size(); i++) + scan[i] = trial(scan_centre.first, scan_centre.second, grid[i], scan_frames); + std::vector peaks; + for (size_t i = 0; i < grid.size(); i++) + if (scan[i] > 0 && (i == 0 || scan[i] >= scan[i - 1]) && (i + 1 == grid.size() || scan[i] >= scan[i + 1])) + peaks.push_back(i); + std::stable_sort(peaks.begin(), peaks.end(), [&](size_t a, size_t b) { return scan[a] > scan[b]; }); + if (peaks.size() > SCAN_PEAKS) + peaks.resize(SCAN_PEAKS); + std::vector> centres{{x, y}, {file_x, file_y}}; centres.insert(centres.end(), candidates.begin(), candidates.end()); float grid_best = file_distance, grid_x = x, grid_y = y; int64_t grid_score = -1; - for (const auto &[cx, cy] : centres) - for (float f = MIN_FACTOR; f <= MAX_FACTOR; f *= STEP) { - const int64_t score = trial(cx, cy, file_distance * f, strongest); - if (score > grid_score) { - grid_score = score; - grid_best = file_distance * f; - grid_x = cx; - grid_y = cy; + for (const size_t p : peaks) + for (size_t i = (p > 0 ? p - 1 : 0); i <= std::min(grid.size() - 1, p + 1); i++) + for (const auto &[cx, cy] : centres) { + const int64_t score = trial(cx, cy, grid[i], strongest); + if (score > grid_score) { + grid_score = score; + grid_best = grid[i]; + grid_x = cx; + grid_y = cy; + } } - } float best_distance = file_distance, best_x = x, best_y = y; int64_t best = score_here; - for (const float f : {0.98f, 0.99f, 1.0f, 1.01f, 1.02f}) { - const int64_t score = trial(grid_x, grid_y, grid_best * f, frames); - if (score > best) { - best = score; - best_distance = grid_best * f; - best_x = grid_x; - best_y = grid_y; + const auto fine = [&](float about, float cx, float cy) { + for (const float f : {0.98f, 0.99f, 1.0f, 1.01f, 1.02f}) { + const int64_t score = trial(cx, cy, about * f, frames); + if (score > best) { + best = score; + best_distance = about * f; + best_x = cx; + best_y = cy; + } } - } + }; + if (!peaks.empty()) + fine(grid_best, grid_x, grid_y); if (score_here < FILE_KEPT_ABOVE * best) { distance = best_distance; x = best_x; @@ -637,8 +781,22 @@ bool Rugnux::DecideStillsGeometry(int start_image, int images_to_process) { experiment_.SetUnitCell(*stills_user_cell_); } logger.Info("Stills detector distance: the file's {:.2f} mm puts {} sampled spots on the lattice, the best " - "between {:.0f}% and {:.0f}% of it {} at {:.2f} mm - the run uses {:.2f} mm", - file_distance, score_here, 100 * MIN_FACTOR, 100 * MAX_FACTOR, best, best_distance, distance); + "between {:.0f} and {:.0f} mm {} at {:.2f} mm - the run uses {:.2f} mm", + file_distance, score_here, MIN_DISTANCE_MM, MAX_DISTANCE_MM, best, best_distance, distance); + // The sample's spots were found at the file's distance, with the spot finder's d-dependent thresholds + // there; from a distance this far off they are found again where the run now is, and the fine step + // asked again on them (a sample found eight times too far out put a seventh as many spots on the + // lattice at the right distance as one found nearer). + if (std::abs(distance / file_distance - 1.0f) > RESAMPLE_ABOVE) { + experiment_.BeamX_pxl(x).BeamY_pxl(y).DetectorDistance_mm(distance); + find_sample(); + tried.clear(); + best = -1; + fine(distance, x, y); + logger.Info("Stills detector distance: sampled again at {:.2f} mm, {} spots on the lattice at {:.2f} mm", + distance, best, best_distance); + distance = best_distance; + } // A centre decided at a distance this far off had little indexing to go on; decided again here. if (distance != file_distance) decide_centre(distance); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 565872458..e7ef757c4 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -1041,6 +1041,8 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { // Geometry overrides (default: keep the value stored in the input file) std::optional beam_x, beam_y, detector_distance_mm, wavelength_A, rot1_rad, rot2_rad, rot3_rad, polarization_factor; + // The input states no detector distance and none was given: the stills search decides it (with -C). + bool detector_distance_unknown = false; std::optional pixel_size_um, gain_per_photon, saturation_value; std::vector extra_masks; // --mask-good / --mask-bad bool detector_mirror_y = false; @@ -1993,8 +1995,10 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { } else if (input_is_xtc) { xtc_reader.ReadFiles(input_files); std::vector missing; - if (!detector_distance_mm) + // Without a distance the stills search finds it, which needs the reference cell (-C). + if (!detector_distance_mm && !fixed_reference_unit_cell) missing.emplace_back("--detector-distance"); + detector_distance_unknown = !detector_distance_mm; if (!xtc_reader.HasPixelSize() && !pixel_size_um) missing.emplace_back("--pixel-size"); if (!xtc_reader.HasWavelength() && !wavelength_A) @@ -2004,7 +2008,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { for (const auto &m: missing) options += (options.empty() ? "" : ", ") + m; logger.Error("The XTC does not state the detector distance{}{}: give {} (the beam centre is " - "measured)", xtc_reader.HasPixelSize() ? "" : ", the pixel size", + "measured; with -C the distance is too)", xtc_reader.HasPixelSize() ? "" : ", the pixel size", xtc_reader.HasWavelength() ? "" : ", the wavelength", options); exit(EXIT_FAILURE); } @@ -3472,6 +3476,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { config.detect_beam_stop = detect_beam_stop; config.beam_center_check = beam_center_check; config.detector_distance_given = detector_distance_mm.has_value(); + config.detector_distance_unknown = detector_distance_unknown; config.detect_saturation = !saturation_value && (input_is_cbf || input_is_marccd || input_is_smv || input_is_crystfel || input_is_panel_table || input_is_cxi || input_is_xtc);