diff --git a/image_analysis/geom_refinement/AssignSpotsToRings.cpp b/image_analysis/geom_refinement/AssignSpotsToRings.cpp index d6152385c..2bc6a0c87 100644 --- a/image_analysis/geom_refinement/AssignSpotsToRings.cpp +++ b/image_analysis/geom_refinement/AssignSpotsToRings.cpp @@ -4,23 +4,42 @@ #include "../../common/JFJochMath.h" #include "AssignSpotsToRings.h" +#include +#include #include #include #include #include -#include #include "../../common/CrystalLattice.h" -FindCircleCenterResult FindCircleCenter(const std::vector &v, int64_t width, int64_t height, int64_t max_spots) { - if ((width <= 0) || (height <= 0) || (max_spots <= 0)) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Invalid image size"); - - std::vector vote(width * height, 0); +FindCircleCenterResult FindCircleCenter(const std::vector &v, int64_t max_spots) { + if (max_spots <= 0) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Invalid spot limit"); + if (v.size() < 3) + return {.total_votes = 0, .votes_for_beam_center = 0, .x = 0.0f, .y = 0.0f}; // Limit to only first 250 spots, given the algorithm is N^3 auto task_size = std::min(v.size(), max_spots); + // The vote grid, one bin per pixel over the spots' bounding box (see the header). + int64_t x0 = std::numeric_limits::max(), x1 = std::numeric_limits::min(); + int64_t y0 = x0, y1 = x1; + for (size_t i = 0; i < task_size; i++) { + x0 = std::min(x0, std::floor(v[i].x)); + x1 = std::max(x1, std::ceil(v[i].x)); + y0 = std::min(y0, std::floor(v[i].y)); + y1 = std::max(y1, std::ceil(v[i].y)); + } + const int64_t width = x1 - x0 + 1; + const int64_t height = y1 - y0 + 1; + if ((width <= 0) || (height <= 0)) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Invalid image size"); + + // uint32_t, not int64_t: the most votes any bin can take is C(max_spots,3), which for the 500-spot + // limit is 2.1e7 - three orders below what a uint32 holds - and it halves the allocation. + std::vector vote(width * height, 0); + for (int i = 0; i < task_size; i++) { for (int j = i+1; j < task_size; j++) { for (int k = j+1; k < task_size; k++) { @@ -36,8 +55,8 @@ FindCircleCenterResult FindCircleCenter(const std::vector &v, int64_ float bx = x1_sq * (v[j].y - v[k].y) + x2_sq * (v[k].y - v[i].y) + x3_sq * (v[i].y - v[j].y); float by = x1_sq * (v[k].x - v[j].x) + x2_sq * (v[i].x - v[k].x) + x3_sq * (v[j].x - v[i].x); - int64_t cx = std::lround(bx / (2.0f * a)); - int64_t cy = std::lround(by / (2.0f * a)); + const int64_t cx = std::lround(bx / (2.0f * a)) - x0; + const int64_t cy = std::lround(by / (2.0f * a)) - y0; if ((cx >= 0) && (cx < width) && (cy >= 0) && (cy < height)) vote[cx + cy * width]++; @@ -63,8 +82,8 @@ FindCircleCenterResult FindCircleCenter(const std::vector &v, int64_ return { .total_votes = total_votes, .votes_for_beam_center = max_votes, - .x = static_cast(cx), - .y = static_cast(cy) + .x = static_cast(cx + x0), + .y = static_cast(cy + y0) }; } @@ -250,27 +269,22 @@ void GuessGeometry(DiffractionGeometry &geom, const std::vector &v, int ring_idx = 0; int cluster_idx = 0; + // Walk the observed clusters and the calibrant's rings together, in ascending q. Which of them + // paired up used to be printed to stdout from here; it is diagnostics a caller wants back as data, + // not something a library should write to a terminal. while (cluster_idx < cluster_annot.size() && ring_idx < ring_q.size()) { - float obs_q = 2 * PI / geom.PxlToRes(cluster_annot[cluster_idx].R_obs); - if (std::fabs(ring_q[ring_idx] - obs_q) < 0.1) { - std::cout << "Found ring " << ring_idx - << " with q_ring = " << ring_q[ring_idx] - << " and q_obs = " << obs_q - << " diff = " << std::fabs(ring_q[ring_idx] - obs_q) << std::endl; + const float obs_q = 2 * PI / geom.PxlToRes(cluster_annot[cluster_idx].R_obs); + if (std::fabs(ring_q[ring_idx] - obs_q) + < RingMatchWindow(ring_q, ring_idx, RING_MATCH_Q_RECIPA)) { for (const auto &spot: cluster_annot[cluster_idx].spots) optimizer_input.push_back({v[spot].x, v[spot].y, ring_q[ring_idx]}); ring_idx++; cluster_idx++; + } else if (ring_q[ring_idx] < obs_q) { + ring_idx++; } else { - std::cout << "Cannot match " << ring_idx - << " with q_ring = " << ring_q[ring_idx] - << " and q_obs = " << obs_q - << " diff = " << std::fabs(ring_q[ring_idx] - obs_q) << std::endl; - if (ring_q[ring_idx] < obs_q) - ring_idx++; - else - cluster_idx++; + cluster_idx++; } } @@ -278,19 +292,33 @@ void GuessGeometry(DiffractionGeometry &geom, const std::vector &v, geom = optimizer.Run(optimizer_input); } +float RingMatchWindow(const std::vector &ring_q, size_t i, float max_window) { + float window = max_window; + if (i > 0) + window = std::min(window, 0.5f * (ring_q[i] - ring_q[i - 1])); + if (i + 1 < ring_q.size()) + window = std::min(window, 0.5f * (ring_q[i + 1] - ring_q[i])); + return window; +} + std::vector AssignSpotsToRings(const DiffractionGeometry &geom, const std::vector &v, const std::vector &ring_q) { std::vector optimizer_input; for (const auto& s: v) { - float q_obs = 2 * PI / geom.PxlToRes(s.x, s.y); - for (const auto &q : ring_q) { - if (std::fabs(q - q_obs) < 0.1) { - optimizer_input.push_back({s.x, s.y, q}); - break; - } + const float q_obs = 2 * PI / geom.PxlToRes(s.x, s.y); + // The NEAREST ring, then the window - not the first ring within a fixed window. Where two rings + // are closer together than that window, taking the first match assigns both of them to the + // lower-q one and biases the distance it fits. + size_t nearest = 0; + for (size_t i = 1; i < ring_q.size(); i++) { + if (std::fabs(ring_q[i] - q_obs) < std::fabs(ring_q[nearest] - q_obs)) + nearest = i; } + if (!ring_q.empty() + && std::fabs(ring_q[nearest] - q_obs) < RingMatchWindow(ring_q, nearest, RING_MATCH_Q_RECIPA)) + optimizer_input.push_back({s.x, s.y, ring_q[nearest]}); } return optimizer_input; } diff --git a/image_analysis/geom_refinement/AssignSpotsToRings.h b/image_analysis/geom_refinement/AssignSpotsToRings.h index 4a806107e..69cd79ef0 100644 --- a/image_analysis/geom_refinement/AssignSpotsToRings.h +++ b/image_analysis/geom_refinement/AssignSpotsToRings.h @@ -21,9 +21,12 @@ struct FindCircleCenterResult { float y; }; +// Beam centre from the circumcentres of spot triples, voted into a pixel grid. The grid spans the +// spots' own bounding box: a powder ring encloses the beam centre, so that is where the answer has to +// be, and it is the only bound the spots themselves justify. It used to be a fixed 4000x4000 box, +// which both allocated 128 MB whatever the detector and put the centre outside the grid - so every +// vote was discarded - on any detector larger than that in either direction. FindCircleCenterResult FindCircleCenter(const std::vector &v, - int64_t width = 4000, - int64_t height = 4000, int64_t max_spots = 500); std::vector> ClusterSpotsIntoRings(const std::vector& r, float eps = 1.5, int minPts = 5); @@ -56,8 +59,19 @@ void GuessGeometry(DiffractionGeometry &geom, const std::vector &v, void OptimizeGeometry(DiffractionGeometry &geom, const std::vector &v, const std::vector &ring_q, bool refine_tilt = true, RingFitUncertainty *unc = nullptr); +// Widest a ring's match window is allowed to be, in q. The clamp below narrows it where the rings +// crowd; this is the value for a well-separated ring. +constexpr float RING_MATCH_Q_RECIPA = 0.1f; + +// Half-width of the window ring i may be matched in, clamped so it can never reach the neighbouring +// ring. Silver behenate's orders sit 0.108 1/A apart and hexagonal ice has three rings within 0.06, so +// a fixed window takes in the neighbour's flank - and a search that stops at its first match then puts +// every point on the lower-q ring of the pair. +float RingMatchWindow(const std::vector &ring_q, size_t i, float max_window); + // Each spot paired with the calibrant ring nearest its observed q, as the points RingOptimizer fits. -// Spots more than 0.1 1/A from every ring are dropped rather than forced onto the closest one. +// Spots further from every ring than RingMatchWindow allows are dropped rather than forced onto the +// closest one. std::vector AssignSpotsToRings(const DiffractionGeometry &geom, const std::vector &v, const std::vector &ring_q); diff --git a/image_analysis/geom_refinement/RingsFromProfile.cpp b/image_analysis/geom_refinement/RingsFromProfile.cpp index b33f26d1c..9402c357c 100644 --- a/image_analysis/geom_refinement/RingsFromProfile.cpp +++ b/image_analysis/geom_refinement/RingsFromProfile.cpp @@ -5,6 +5,7 @@ #include #include "RingsFromProfile.h" +#include "AssignSpotsToRings.h" // RingMatchWindow #include "../../common/JFJochMath.h" namespace { @@ -22,6 +23,15 @@ float SectorPeakQ(const std::vector &profile, int32_t q_bins, int phi_bin const auto value = [&](int i) { return profile[row + static_cast(i)]; }; const auto q_of = [&](int i) { return low_q + (static_cast(i) + 0.5f) * q_spacing; }; + // A bin no pixel fell in is NaN, not zero (AzimuthalIntegrationProfile::GetResult), and the four + // background bins are where a module gap or the beam stop shows up first. Say so rather than + // relying on NaN comparisons to fail the peak test further down: a sector whose background cannot + // be measured has no measurable peak either. + for (int i : {lo_bin, lo_bin + 1, hi_bin - 1, hi_bin}) { + if (!std::isfinite(value(i))) + return NAN; + } + const float bkg_lo = 0.5f * (value(lo_bin) + value(lo_bin + 1)); const float bkg_hi = 0.5f * (value(hi_bin) + value(hi_bin - 1)); const auto bkg_at = [&](int i) { @@ -99,11 +109,7 @@ std::vector RingsFromAzimuthalProfile(const std::vector 0) - window = std::min(window, 0.5f * (q_ring - calibrant_ring_q[i - 1])); - if (i + 1 < calibrant_ring_q.size()) - window = std::min(window, 0.5f * (calibrant_ring_q[i + 1] - q_ring)); + const float window = RingMatchWindow(calibrant_ring_q, i, q_window_recipA); if (!(q_ring - window > low_q) || !(q_ring + window < high_q)) continue; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 3e3c661be..dd3465172 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -1988,6 +1988,12 @@ static int RunRugnux(int argc, char **argv) { config.fit_spindle = fit_spindle; config.write_process_h5 = false; // the .poni below is the output of this mode + // Ice-ring handling has no place in a calibration run and does active harm. Flagged spots are + // sorted LAST by the spot budget (FilterSpotsByCount), so they are the first thrown away - which + // for --calibrant ice discards exactly the spots being calibrated on, and for every other + // standard discards whichever of its rings happen to fall in the fixed hexagonal bands. + experiment.DetectIceRings(false); + // Spot finding for --calibration spots. Indexing is off: a calibration wants the spot positions // and nothing else, and the calibrant is a powder with no lattice to index. config.spot_finding.enable = true; diff --git a/tests/CalibrationTest.cpp b/tests/CalibrationTest.cpp index 2226b7a37..1f6bda6c4 100644 --- a/tests/CalibrationTest.cpp +++ b/tests/CalibrationTest.cpp @@ -10,6 +10,7 @@ #include "../common/Definitions.h" #include "../common/JFJochMath.h" +#include "../image_analysis/geom_refinement/AssignSpotsToRings.h" #include "../image_analysis/geom_refinement/Calibrants.h" #include "../image_analysis/geom_refinement/PowderCalibration.h" @@ -128,3 +129,42 @@ TEST_CASE("Calibration_PoniFileAxisConvention", "[DetGeomCalib]") { + std::to_string(x.GetXPixelsNumConv()) + "]"; CHECK(keys["Detector_config"].find(shape) != std::string::npos); } + +// The match window may never reach the neighbouring ring, for any calibrant. Where two rings are closer +// together than twice the nominal window, a fixed window takes in the neighbour's flank - which the +// rings path reads as this ring's background, and which the spots path (before it took the NEAREST ring) +// resolved by assigning both to the lower-q one. +TEST_CASE("Calibration_RingMatchWindowNeverReachesTheNeighbour", "[DetGeomCalib]") { + for (const auto &c : Calibrants()) { + const auto q = CalibrantRings(c.name); + REQUIRE(q.size() > 1); + for (size_t i = 0; i < q.size(); ++i) { + const float w = RingMatchWindow(q, i, RING_MATCH_Q_RECIPA); + CHECK(w <= RING_MATCH_Q_RECIPA); + CHECK(w > 0.0f); + if (i > 0) + CHECK(q[i] - w >= 0.5f * (q[i] + q[i - 1]) - 1e-6f); + if (i + 1 < q.size()) + CHECK(q[i] + w <= 0.5f * (q[i] + q[i + 1]) + 1e-6f); + } + } +} + +// ...and the clamp is not a no-op. Silver behenate's orders sit about 0.108 1/A apart and hexagonal ice +// has rings inside 0.06, so both are narrowed below the nominal window - while LaB6, whose rings are +// well separated at low q, keeps it. Without a standard that actually crowds, the test above would pass +// on a clamp that never fired. +TEST_CASE("Calibration_CrowdedStandardsNarrowTheWindow", "[DetGeomCalib]") { + auto narrowed = [](const std::string &name) { + const auto q = CalibrantRings(name); + size_t n = 0; + for (size_t i = 0; i < q.size(); ++i) + if (RingMatchWindow(q, i, RING_MATCH_Q_RECIPA) < RING_MATCH_Q_RECIPA) ++n; + return n; + }; + CHECK(narrowed("agbh") > 0); + CHECK(narrowed("ice") > 0); + // The innermost LaB6 rings are more than 0.2 1/A apart, so nothing narrows them. + const auto lab6 = CalibrantRings("lab6"); + CHECK(RingMatchWindow(lab6, 0, RING_MATCH_Q_RECIPA) == Catch::Approx(RING_MATCH_Q_RECIPA)); +}