diff --git a/image_analysis/bragg_integration/SpotFootprint.cpp b/image_analysis/bragg_integration/SpotFootprint.cpp index bfa6d0257..584d973da 100644 --- a/image_analysis/bragg_integration/SpotFootprint.cpp +++ b/image_analysis/bragg_integration/SpotFootprint.cpp @@ -5,6 +5,8 @@ #include #include +#include +#include namespace { @@ -124,3 +126,50 @@ SpotFootprint FootprintFromSpots(const std::vector &spots, float } return fp; } + +void MeasureFootprintOffsets(const std::vector &spots, const std::vector &reflections, + float beam_x, float beam_y, std::vector &out) { + std::map, const Reflection *> predicted; + for (const auto &r : reflections) + predicted[{r.h, r.k, r.l}] = &r; + for (const auto &s : spots) { + if (!s.indexed || s.lattice != 0) continue; + const auto it = predicted.find({static_cast(s.h), static_cast(s.k), static_cast(s.l)}); + if (it == predicted.end()) continue; + const float px = it->second->predicted_x, py = it->second->predicted_y; + const float rx = px - beam_x, ry = py - beam_y; + const float r = std::sqrt(rx * rx + ry * ry); + if (!(r > 1.0f)) continue; + const float ux = rx / r, uy = ry / r; + const float dx = s.x - px, dy = s.y - py; + out.push_back({r, dx * ux + dy * uy, -dx * uy + dy * ux}); + } +} + +SpotFootprint FootprintWithOffsets(const SpotFootprint &widths, std::vector offsets) { + if (widths.empty()) return widths; + // Sorted, so the sums below do not depend on the order the frames were measured in. + std::sort(offsets.begin(), offsets.end()); + const int n = static_cast(widths.sigma_rad.size()); + std::vector s2r(n, 0.0), s2t(n, 0.0); + std::vector cnt(n, 0); + for (const auto &o : offsets) { + const int b = std::clamp(static_cast(o.r_px / widths.bin_px), 0, n - 1); + s2r[b] += static_cast(o.off_rad) * o.off_rad; + s2t[b] += static_cast(o.off_tan) * o.off_tan; + ++cnt[b]; + } + std::vector filled; + for (int b = 0; b < n; ++b) + if (cnt[b] >= FOOTPRINT_MIN_SPOTS_PER_BIN) filled.push_back(b); + if (filled.empty()) return widths; + SpotFootprint fp = widths; + for (int b = 0; b < n; ++b) { + int best = filled.front(); + for (int f : filled) + if (std::abs(f - b) < std::abs(best - b)) best = f; + fp.sigma_rad[b] = static_cast(std::sqrt(widths.sigma_rad[b] * widths.sigma_rad[b] + s2r[best] / cnt[best])); + fp.sigma_tan[b] = static_cast(std::sqrt(widths.sigma_tan[b] * widths.sigma_tan[b] + s2t[best] / cnt[best])); + } + return fp; +} diff --git a/image_analysis/bragg_integration/SpotFootprint.h b/image_analysis/bragg_integration/SpotFootprint.h index d973481a2..9b54c8186 100644 --- a/image_analysis/bragg_integration/SpotFootprint.h +++ b/image_analysis/bragg_integration/SpotFootprint.h @@ -23,12 +23,25 @@ // to the integrator through BraggIntegrationSettings; where it says a spot outgrows the r1 disk, the // background ring is moved clear of the spot and the profile is fitted at the measured width // (BraggStencil.h). Compact spots leave the integration exactly as it was. +// +// A spot can also be in the wrong PLACE for the disk: a crystal of slightly misaligned domains - the +// ferroelastic domains of a crystal below a phase transition, a cracked or split crystal - records +// each reflection as two or more spots around the position the averaged lattice predicts, and they +// move apart with resolution. Each of them is compact, so the widths above, measured about each spot, +// say nothing; the r1 disk holds the gap between them and the background ring lands on them. What +// does see it is where the spots sit against the prediction: once the lattice is known, every +// indexed spot is compared with the predicted position of its own reflection on the same frame, and +// the mean square of that offset is added to the width of the spot, radially and tangentially apart. +// The table then describes the reflection rather than the spot - and on a crystal whose spots sit on +// their predictions it moves the widths by the prediction error alone, a fraction of a pixel. // ============================================================================= #include #include #include "../../common/BraggIntegrationSettings.h" +#include "../../common/Reflection.h" +#include "../../common/SpotToSave.h" // One spot's measured widths, and its distance from the beam centre [px]. struct FootprintSpot { @@ -48,6 +61,27 @@ void MeasureFootprintSpots(const int32_t *img, int width, int height, float beam // holding too few spots takes the nearest bin that has enough; with no such bin the table is empty. SpotFootprint FootprintFromSpots(const std::vector &spots, float r_max); +// One indexed spot's offset from the predicted position of its reflection [px], along and across the +// radius, and the prediction's distance from the beam centre [px]. +struct FootprintOffset { + float r_px = 0.0f; + float off_rad = 0.0f; + float off_tan = 0.0f; + bool operator<(const FootprintOffset &o) const { + return r_px != o.r_px ? r_px < o.r_px : off_rad != o.off_rad ? off_rad < o.off_rad : off_tan < o.off_tan; + } +}; + +// The offsets of one frame: every spot indexed on the main lattice against the reflection of the same +// hkl predicted on that frame. Spots whose reflection is not among the predictions are left out. +void MeasureFootprintOffsets(const std::vector &spots, const std::vector &reflections, + float beam_x, float beam_y, std::vector &out); + +// The width table with the offsets added in quadrature: in each distance bin of `widths`, the mean +// square offset along and across the radius. A bin holding too few offsets takes the nearest bin +// that has enough; with no such bin, or no widths, the widths come back unchanged. +SpotFootprint FootprintWithOffsets(const SpotFootprint &widths, std::vector offsets); + // Spots per distance bin a table needs before it is believed. constexpr int FOOTPRINT_MIN_SPOTS_PER_BIN = 20; constexpr int FOOTPRINT_BINS = 12; diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 596ee0b0e..d8c93cde8 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -1311,6 +1311,8 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru const float H = static_cast(experiment_.GetYPixelsNum()); const SpotFootprint fp = FootprintFromSpots(footprint_spots, std::hypot(std::max(bx, W - bx), std::max(by, H - by))); + if (!fp.empty()) + prescan_footprint_ = fp; // It acts only where a spot outgrows the r1 disk (BraggStencil.h), so a pattern of compact // spots is left on exactly the settings - and the passes - it had without it. const float r1_now = experiment_.GetBraggIntegrationSettings().GetR1(); @@ -3719,8 +3721,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b indexing_bound_A = indexing_settings.GetFFT_MaxUnitCell_A(); indexer_pool = std::make_unique(indexing_settings, IndexerConstruction::OnFirstUse); indexer = std::make_unique(experiment_, indexer_pool.get()); - // With no per-image file to write, nothing reads the message's reflection list. - if (!writer_queue) + // With no per-image file to write, nothing reads the message's reflection list - except the + // geometry pre-pass, which measures the spots against it (SpotFootprint.h). + if (!writer_queue && !(geometry_prepass && prescan_footprint_)) indexer->KeepReflectionsInMessage(false); // The reference that breaks the per-image indexing ambiguity: a reference MTZ where there is // one, otherwise intensities computed from the model, which carry the same information. Serial @@ -6117,6 +6120,12 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b std::mutex harmonic_m; HarmonicEvidence harmonic_evidence; + // Where the indexed spots sit against their predicted reflections, measured on the geometry + // pre-pass for the canonical pass's footprint (SpotFootprint.h). + const bool measure_offsets = geometry_prepass && prescan_footprint_.has_value(); + std::mutex footprint_m; + std::vector footprint_offsets; + auto azint_worker = [&]() { std::vector decompression_buffer; std::shared_ptr img; @@ -6192,6 +6201,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b /*enable_fused_adaptive_gpu=*/true); AzimuthalIntegrationProfile profile(mapping); HarmonicEvidence worker_harmonic; + std::vector worker_offsets; while (!cancelled_) { const int ordinal = next_ordinal.fetch_add(1); @@ -6240,6 +6250,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b } AddHarmonicEvidence(msg.spots, experiment_.GetDiffractionGeometry(), worker_harmonic); + if (measure_offsets) + MeasureFootprintOffsets(msg.spots, msg.reflections, experiment_.GetBeamX_pxl(), + experiment_.GetBeamY_pxl(), worker_offsets); plots.Add(msg, profile); if (writer_queue) writer_queue->Post(msg, img); @@ -6253,6 +6266,10 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b std::lock_guard lock(harmonic_m); harmonic_evidence.Add(worker_harmonic); } + { + std::lock_guard lock(footprint_m); + footprint_offsets.insert(footprint_offsets.end(), worker_offsets.begin(), worker_offsets.end()); + } const auto bragg = analysis.BraggCounts(); bragg_predicted += bragg.predicted; @@ -6306,6 +6323,31 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b result.cancelled = cancelled_; result.images_processed = finished_count.load(); + // The footprint of the reflection rather than of the spot: the pre-scan's widths plus where the + // spots sit against their predictions (SpotFootprint.h). Like the widths it acts only where the + // reflection outgrows the r1 disk, and like them it is held back for the canonical pass. + if (measure_offsets && !cancelled_) { + const SpotFootprint fp = FootprintWithOffsets(*prescan_footprint_, footprint_offsets); + const BraggIntegrationSettings canonical = bragg_adaptive_.value_or(experiment_.GetBraggIntegrationSettings()); + bool outgrows = false; + for (size_t b = 0; b < fp.sigma_rad.size(); ++b) + outgrows |= BRAGG_FOOTPRINT_NSIGMA * std::max(fp.sigma_rad[b], fp.sigma_tan[b]) > canonical.GetR1(); + std::string table; + for (size_t b = 0; b < fp.sigma_rad.size(); ++b) + table += fmt::format(" {:.0f}:{:.1f}/{:.1f}", (b + 0.5f) * fp.bin_px, fp.sigma_rad[b], fp.sigma_tan[b]); + logger.Info("Reflection footprint: the spot widths plus the offsets of {} indexed spots from their " + "predicted reflections, sigma radial/tangential [px] by distance from the beam [px]:{}{}", + footprint_offsets.size(), table, + outgrows ? "" : fmt::format(" - every reflection fits the r1={:.1f} disk", canonical.GetR1())); + if (outgrows) { + if (!bragg_before_adaptive_) + bragg_before_adaptive_ = experiment_.GetBraggIntegrationSettings(); + BraggIntegrationSettings bis = canonical; + bis.Footprint(fp); + bragg_adaptive_ = bis; + } + } + if (full && indexer && experiment_.IsRotationIndexing() && !cancelled_) { const SupercellProbe probe = indexer->GetSupercellProbe(); if (probe[0][0].n > 0 && probe[0][1].n > 0) { diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 638b7502c..e4be283b1 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -719,6 +719,11 @@ class Rugnux { // radius before the pre-pass and applying it after keeps the geometry - and so the lattice - exactly // what the fixed radius produced, and still integrates the canonical pass at the measured width. std::optional bragg_adaptive_; + // The spot footprint as the pre-scan measured it, about each spot (SpotFootprint.h), whether or + // not it was adopted. The geometry pre-pass adds to it how far the spots sit from their predicted + // reflections - which takes a lattice the pre-scan does not have - and the sum goes to the + // canonical pass in bragg_adaptive_. Kept apart so a repeated pre-pass starts from it again. + std::optional prescan_footprint_; // Fraction of the predicted reflections the last completed pass dropped because NEIGHBOURING // reflections left their background ring below six clean pixels. Measured, not predicted - it is // the only thing that tells a radius the pattern can take from one it cannot (see RunAllPasses). diff --git a/tests/SpotFootprintTest.cpp b/tests/SpotFootprintTest.cpp index 46699b50e..c74d65e2d 100644 --- a/tests/SpotFootprintTest.cpp +++ b/tests/SpotFootprintTest.cpp @@ -66,3 +66,39 @@ TEST_CASE("SpotFootprint_TableFillsSparseBinsFromNeighbours", "[Integration][por REQUIRE(fp.sigma_rad[FOOTPRINT_BINS - 3] == 3.0f); // nearer the last REQUIRE(FootprintFromSpots({}, 1200.0f).empty()); } + +// A reflection recorded as a doublet: two indexed spots on either side of the prediction, along the +// radius. The offsets are found against the reflection of the same hkl, and their mean square is +// added to the widths of the bin they fall in. +TEST_CASE("SpotFootprint_OffsetsFromPredictionWidenTheTable", "[Integration][portable]") { + const float bx = 600.0f, by = 600.0f; + std::vector refl(1); + refl[0].h = 1; refl[0].k = 2; refl[0].l = 3; + refl[0].predicted_x = bx + 1000.0f; // on +x, so radial = x + refl[0].predicted_y = by; + std::vector spots; + spots.push_back({.x = bx + 1004.0f, .y = by, .lattice = 0, .h = 1, .k = 2, .l = 3, .indexed = true}); + spots.push_back({.x = bx + 996.0f, .y = by, .lattice = 0, .h = 1, .k = 2, .l = 3, .indexed = true}); + spots.push_back({.x = bx + 1000.0f, .y = by + 30.0f, .lattice = 0, .h = 3, .k = 2, .l = 1, .indexed = true}); // no such prediction + spots.push_back({.x = bx + 1000.0f, .y = by + 30.0f, .lattice = -1, .h = 1, .k = 2, .l = 3, .indexed = false}); // not indexed + std::vector offsets; + MeasureFootprintOffsets(spots, refl, bx, by, offsets); + REQUIRE(offsets.size() == 2); + REQUIRE_THAT(offsets[0].off_rad, Catch::Matchers::WithinAbs(4.0, 1e-4)); + REQUIRE_THAT(offsets[0].off_tan, Catch::Matchers::WithinAbs(0.0, 1e-4)); + + std::vector pool; + for (int i = 0; i < FOOTPRINT_MIN_SPOTS_PER_BIN; ++i) + pool.insert(pool.end(), offsets.begin(), offsets.end()); + SpotFootprint widths; + widths.bin_px = 100.0f; + widths.sigma_rad.assign(12, 1.0f); + widths.sigma_tan.assign(12, 2.0f); + const SpotFootprint fp = FootprintWithOffsets(widths, pool); + for (int b = 0; b < 12; ++b) { // one filled bin (10) serves them all + REQUIRE_THAT(fp.sigma_rad[b], Catch::Matchers::WithinAbs(std::sqrt(17.0), 1e-4)); + REQUIRE_THAT(fp.sigma_tan[b], Catch::Matchers::WithinAbs(2.0, 1e-4)); + } + REQUIRE(FootprintWithOffsets(widths, {}).sigma_rad[3] == 1.0f); + REQUIRE(FootprintWithOffsets(SpotFootprint{}, pool).empty()); +}